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ABSTRACT 


Transitional millisecond pulsars are an emerging class of sources linking low-mass X-ray binaries to millisecond radio pulsars in 
binary systems. These pulsars alternate between a radio pulsar state and an active low-luminosity X-ray disc state. During the active 
state, these sources exhibit two distinct emission modes (high and low) that alternate unpredictably, abruptly, and incessantly. X-ray 
to optical pulsations are observed only during the high mode. Knowledge of the root reason for this puzzling behaviour remains 
elusive. This paper presents the results of the most extensive multi-wavelength campaign ever conducted on the transitional pulsar 
prototype, PSR J1023--0038, covering from radio to X-rays. The campaign was carried out over two nights in June 2021, and involved 
12 different telescopes and instruments including XMM-Newton, HST, VLT/FORS2 (in polarimetric mode), ALMA, VLA and FAST. 
By modelling the broadband spectral energy distributions in both emission modes, we show that the mode switches are caused by 
changes in the innermost region of the accretion disc. These changes trigger the emission of discrete mass ejections, which occur on 
top of a compact jet, as testified by the detection of at least one short-duration millimetre flare with ALMA at the high-to-low mode 
switch. A subsequent re-enshrouding of the pulsar completes the scenario behind the mode switches. 


Key words. accretion, accretion disks — ISM: jets and outflows — pulsars: general — radio continuum: stars — stars: neutron — X-rays: 
binaries 


1. Introduction 


The evolutionary path that leads a neutron star (NS) in a low- 


7 : . 
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passes through different stages, the most enigmatic of which is 
represented by transitional millisecond pulsars (tMSPs). These 
sources alternate between a radio pulsar state and an active 
low-luminosity X-ray disc state on timescales from weeks to 
months (for reviews, see Campana & Di Salvo 2018; Papitto 
& de Martino 2022). The archetypal tMSP is PSR J1023+0038 
(hereafter referred to as J1023), which was discovered in 2007 
as a 1.69 ms radio pulsar orbiting a low-mass (20.2 Mo) com- 
panion star with a period of 4.75 hr (Archibald et al. 2009). In 
2013, J1023 showed a sudden brightening in the emission at X- 
ray and gamma-ray frequencies (by a factor of 5—10) as well as 
at ultraviolet (UV) and optical frequencies (by 1—2 mag), along 
with the disappearance of the pulsed radio signal (Stappers et al. 
2014; Patruno et al. 2014). Double-peaked optical emission lines 
were soon detected, indicating the formation of an accretion disc 
(Coti Zelati et al. 2014). J1023 has remained in this active state 
ever since, at an X-ray luminosity of Ly = 7 x 10? erg s^! in the 
0.3-80 keV energy range (Coti Zelati et al. 2018) for a distance 
of 1.37 kpc (Deller et al. 2012). 

During the sub-luminous X-ray state, the X-ray emission 
of J1023 is observed to switch between two distinct intensity 
modes, dubbed high and low, with sporadic flares in between. 
The high mode occurs 70-80% of the time, while the low mode 
occurs 20-30% of the time and typically lasts from a few tens 
of seconds to minutes. The drop and rise timescales involved in 
the mode switching are of the order of «10 s. X-ray, UV and op- 
tical pulsations occur simultaneously in the X-ray high-intensity 
mode, and disappear in the low mode (Archibald et al. 2015; Pa- 
pitto et al. 2019; Miraval Zanon et al. 2022; Illiano et al. 2023). 
The UV, optical, and near infrared emission is dominated by the 
accretion disc and irradiated companion star and also shows flar- 
ing and flickering activities (Shahbaz et al. 2015; Kennedy et al. 
2018; Papitto et al. 2018; Hakala & Kajava 2018; Shahbaz et al. 
2018). Optical polarimetric observations reveal a linear polarisa- 
tion at a level of ~1% (Baglio et al. 2016; Hakala & Kajava 2018) 
of uncertain origin. Bright, variable radio continuum emission is 
also detected (Deller et al. 2015). In particular, the low X-ray 
mode episodes are accompanied by an increase of the radio flux 
by a factor of «3 on average, with fairly equal rise and decay 
times and an overall duration that matches low X-ray mode in- 
tervals (Bogdanov et al. 2018, see also Fig. 1). These rebrighten- 
ings are associated with a temporal evolution of the radio spec- 
trum from inverted to steep. Currently, J1023 spins down at a 
rate that is ~4% faster than in the radio pulsar state (Burtovoi 
et al. 2020). This overall behaviour is unique in the landscape of 
X-ray binaries and still defies a comprehensive explanation. 

This paper presents the results of the largest multi- 
wavelength campaign ever performed on J1023. Our aim is 
to ultimately understand the root cause of the mode changes 
in the active sub-luminous X-ray state and the implication for 
our understanding of the accretion-ejection coupling in X-ray 
binary systems harboring rapidly spinning NSs. The observa- 
tions of our campaign were conducted on two different nights 
in June 2021. The first night (Night 1) involved XMM-Newton, 
the Nuclear Spectroscopic Telescope Array (NuSTAR), the Hub- 
ble Space Telescope (HST), the Son OF ISAAC (SOFI) in- 
strument mounted on the New Technology Telescope (NTT), 
the Karl G. Jansky Very Large Array (VLA) and the Five- 
hundred-meter Aperture Spherical Telescope (FAST); the sec- 
ond night (Night2) involved the Neutron Star Interior Compo- 
sition Explorer Mission (NICER), the Neil Gehrels Swift Ob- 
servatory (Swift), the Very Large Telescope (VLT), the Rapid 
Eye Mount (REM) telescope and the Atacama Large millime- 
tre/submillimetre Array (ALMA). The paper is structured as fol- 
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lows: in Sect.2, we present the observations and data analysis. 
The results are then reported in Sect. 3. In Sect. 4, we propose 
a scenario for the mode switching, present the results of the 
modelling of the spectral energy distributions (SEDs) in the two 
modes, and discuss the results in relation to our scenario. Con- 
clusions follow in Sect. 5. 


2. Observations and data analysis 


Table 1 lists all the observations of J1023 presented in this paper. 
In the following, we describe these observations and the proce- 
dures adopted to process and analyse these data. 


2.1. XMM-Newton (Night 1) 


XMM-Newton (Jansen et al. 2001; Schartel et al. 2022) observed 
J1023 on 2021 June 3-4 (Obs ID. 0864010101) using the Eu- 
ropean Photon Imaging Cameras (EPIC) and the Optical/UV 
Monitor telescope (OM; Mason et al. (2001)). The EPIC-pn 
(Strüder et al. 2001) was set in fast timing mode (time res- 
olution of 29.52 us), the two EPIC-MOS (Turner et al. 2001) 
in full frame mode (time resolution of 2.6s) and the OM in 
the fast window mode (time resolution of 0.5 s) with the B fil- 
ter (effective wavelength 4330 À; full-width at half-maximum 
926 À). The Observation Data Files were processed and anal- 
ysed using the Science Analysis System (SAS; v. 20.0). The X- 
ray data were analysed in the same way as described by Miraval 
Zanon et al. (2022). No periods of high background activity were 
detected. The background-subtracted light curve was extracted 
over the time interval during which all three EPIC instruments 
collected data simultaneously and was binned with a time res- 
olution of 10s. The optical background-subtracted light curve 
was extracted using the oMrcHaiN pipeline with default parame- 
ters. The light curve was binned at a time resolution of 120s to 
ensure a high SNR and detect any changes in the optical bright- 
ness at least during prolonged periods of low X-ray mode. The 
net count rates were then converted to magnitudes in the Vega 
system using the formula provided on the XMM-Newton online 
threads!. To align the X-ray and optical light curves with those 
extracted from other observatories, the arrival times of photons 
were converted from the Terrestrial Time (TT) standard to the 
Coordinated Universal Time (UTC) standard without applying a 
barycentric correction’. 


2.2. NuSTAR (Night 1) 


NuSTAR (Harrison et al. 2013) observed J1023 for a total of 
—48.6 ks and a net exposure of ~22 ks (Obs ID. 30601005002), 
almost fully overlapping with the XMM-Newton observation. 
Data were processed and analysed using the NuSTAR Data 
Analysis Software (NUSTARDAS, v.2.1.2) with the instrumen- 
tal calibration files stored in CALDB v20220525 and the same 
prescriptions reported by Miraval Zanon et al. (2022). J1023 was 
detected up to energies of ~50 keV in both focal plane modules 
(FPMA and FPMB). Light curves were extracted separately for 
the two modules, combined to increase the SNR, and binned at 
a time resolution of 50 s. 


! https://www.cosmos.esa.int/web/xmm-newton/ 
sas-watchout-uvflux 

? This same procedure was also applied to the data collected with Nus- 
TAR, Swift and NICER presented in the following. 
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Table 1. Observation Log of J1023 in June 2021. Observations are listed in order of decreasing observing frequency. 


Telescope/Instrument Obs. ID/Project Start — End time Exposure Band (setup) 
MMM DD hh:mm:ss (UTC) (ks) 

Night 1 
XMM-Newton/EPIC-pn 0864010101 Jun 3 20:44:32 — Jun 4 13:39:15 58.7 0.3-10 keV (fast timing) 
XMM-Newton/EPIC-MOS1 0864010101 Jun 3 21:09:37 — Jun 4 13:39:01 58.1 0.3-10 keV (full frame) 
XMM-Newton/EPIC-MOS2 0864010101 Jun 3 22:31:54 — Jun 4 13:39:02 53.5 0.3-10 keV (full frame) 
NuSTAR/FPMA/FPMB 30601005002 Jun 3 20:36:09 — Jun 4 10:06:09 22.3/22.1 3—79 keV 
Swift/XRT 00033012207 Jun 3 18:29:46 — Jun 3 23:45:44 6.4 0.3-10 keV (photon counting) 
Swift/XRT 00033012208 Jun 4 02:27:57 — Jun 4 12:17:45 9.5 0.3-10 keV (photon counting) 
Swift/XRT 00033012209 Jun 4 03:50:46 — Jun 4 08:44:56 2.0 0.3-10 keV (photon counting) 
SwifiIUVOT 00033012207 Jun 3 18:29:50 — Jun 3 23:45:46 6.3 UVM2 (Event) 
Swift/U VOT 00033012209 Jun 4 03:50:50 — Jun 4 08:44:58 1.9 UVM2 (Event) 
HST STIS/NUV-MAMA 16061 Jun 4 00:09:07 — Jun 4 00:43:27 2.1 G230L (fast timing spectroscopy) 
HST STIS/NUV-MAMA 16061 Jun 4 01:35:28 — Jun 4 02:20:18 2.7 G230L (fast timing spectroscopy) 
XMM-Newton/OM 0864010101 Jun 3 22:51:01 — Jun 4 13:23:52 48.4 B (fast window) 
NTT/SOFI Jun 3 23:49:09 — Jun 4 02:42:38 10.4 J (fast photometry) 
VLA SJ6401 Jun 3 22:30:00 - Jun 4 02:29:20 6.7 Band X (C configuration) 
FAST PT2020_0044 Jun 4 08:47:45 — Jun 4 11:47:45 10.8 Band L 

Night 2 
NICER/XTI 4034060102 Jun 26 21:56:04 — Jun 26 22:34:31 2.0 0.3-10 keV 
NICER/XTI 4034060103 Jun 26 23:29:04 — Jun 27 01:33:51 4.1 0.3-10 keV 
Swift/XRT 00033012214 Jun 26 22:19:35 — Jun 26 22:39:47 1.2 0.3-10 keV (photon counting) 
Swift/XRT 00033012213 Jun 26 23:56:12 — Jun 26 23:58:46 0.15 0.3-10 keV (photon counting) 
Swift/XRT 00033012215 Jun 27 00:04:44 — Jun 27 00:21:47 1.0 0.3-10 keV (photon counting) 
Swift/UVOT 00033012214 Jun 26 22:19:35 — Jun 26 22:39:47 1.2 UVM2 (Event+Imaging) 
Swift/{UVOT 00033012213 Jun 26 23:56:12 — Jun 26 23:58:46 0.15 UVM2 (Event+Imaging) 
SwifiIUVOT 00033012215 Jun 27 00:04:44 — Jun 27 00:21:47 1.0 UVM2 (Event+Imaging) 
VLT/FORS2 107.22RK.001 Jun 26 22:51:56 — Jun 27 01:15:30 8.7 R 
REM/ROSS2 43013 Jun 26 23:20:33 — Jun 27 01:24:41 7.4 griz 
REM/REMIR 43013 Jun 26 23:20:44 — Jun 27 01:25:34 7.4 K 
ALMA 2019.A.00036.S Jun 26 21:58:42 — Jun 27 01:12:06 9.3 Band 3 (C43-7 configuration) 


2.3. HST (Night 1) 


The Space Telescope Imaging Spectrograph (STIS; Woodgate 
et al. 1998) on board HST observed J1023 for 4750s using 
the near UV Multi-Anode Micro-channel Array detector (NUV- 
MAMA) in TIME-TAG mode (time resolution of 125 us); the 
data were collected with the G230L grating equipped with a 52" 
X 0.2 slit with a spectral resolution of ~500 over the nominal 
range (first order). The data analysis procedure that we adopted 
is the same as that described by Miraval Zanon et al. (2022). In 
summary, we used the stis photons package? to correct the 
position of the slit channels and to assign the wavelengths to each 
time of arrival (ToA). We selected ToAs belonging to channels 
993-1005 of the slit and in the 165—310 nm wavelength interval 
to isolate the source signal. 


2.4. NTT/SOFI (Night 1) 


We collected NIR high time-resolution imaging data using the 
SOFI (Moorwood et al. 1998) instrument mounted on the Euro- 
pean Southern Observatory (ESO) 3.6-m NTT at La Silla Obser- 
vatory (Chile). Observations were carried out using the J filter 
from 2021 June 3 at 23:49:09 UTC to June 4 at 02:42:38 UTC 
(Table 1). In order to achieve 1-s time resolution, observations 
were performed in Fast Photometry mode, which restricts the 
area on which the detector is read. Specifically, the instrument 
was set to read a 400x400 pixels window (1.9' x 1.9’). During 
the observation, a series of 100 frames each with an exposure 
time of 1.0 s were stacked together to form a *data-cube". There 
was a gap of ~3s between each cube. Photometric data were 


3 https: //github.com/Alymantara/stis_photons 


extracted through aperture photometry with an aperture radius 
of 1.4" using a custom pipeline based on the DAOPHOT package 
(Stetson 1987b). In order to minimise systematic errors, we per- 
formed photometry on J1023 and a reference star (at position 
R.A. = 10^23750554, Dec. = --00*38'16."4, J2000.0, and with a 
magnitude of J = 15.44 + 0.05 mag in the Vega system from the 
Two Micron All-Sky Survey (2MASS) point source catalogue), 
assumed to have constant brightness. We then re-binned the light 
curve at a time resolution of 10 s to reduce the scatter and investi- 
gate the presence of variability at the same time resolution as the 
light curves in the X-ray and UV bands. For a more conservative 
approach, we evaluated the errors on the magnitudes of the re- 
binned light curves as the standard deviation of the magnitudes 
measured in each bin. 


2.5. VLA (Night 1) 


The VLA observed J1023 for 4 hours, starting on 2021 June 3 
at 22:30 UTC. The array was in its C-configuration and data 
were taken at a central frequency of 10 GHz, with a bandwidth 
of 4 GHz. J1331+3030 was used as flux and bandpass calibra- 
tor, while J1024—0052 was used as phase calibrator. The separa- 
tion between the target and the phase calibrator is 0.3°. To cali- 
brate and analyse the data we followed the prescription of Deller 
et al. (2015) and Bogdanov et al. (2018). Here we report the main 
steps. After initial calibration using the VLA pipeline within the 
Common Astronomy Software Application (cAsA) package (ver- 
sion 5.1.1), manual flagging of bad data was performed. Unfortu- 
nately, the first half of the observation turned out to be unusable, 
possibly due to the bad weather conditions that occurred during 
the observation, and was completely flagged out. 
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Subsequently, we self-calibrated the data using a solution 
interval of 75s for the phase calibration and 45 min for the 
amplitude+phase calibration. Upon inspection of the field of 
view, three main sources were identified: the target J1023, 
a faint source to the north, namely J102348.2+004017, and 
the brightest source in the field, located to the south-east, i.e. 
J102358.2+003826. J102358.2+003826 is a well-known galaxy 
at redshift z — 0.449 (Ahn et al. 2014) with a resolved ex- 
tended emission whose secondary lobes heavily affect the over- 
all rms noise level (Deller et al. 2015). To remove the prob- 
lematic contribution of J102358.2--003826, we first cleaned the 
image using a mask with the tclean task. We set the nterms 
parameter equal to 2 to take into account the spectral proper- 
ties of the galaxy over the wide bandwidth of the observation 
(Deller et al. 2015). This step provided the model visibilities of 
J102358.2+003826. Next, we used the uvsub task to subtract 
these model visibilities from the self-calibrated visibilities and 
deconvolved the residuals. The resulting image therefore does 
not have any contribution from J102358.2+003826. 

As pointed out by Deller et al. (2015), the gain solu- 
tions from the self-calibration obtained before the subtraction 
of J102358.2+003826 were dominated by the pointing errors 
of this source. This source was located beyond the half-power 
point of the antenna, resulting in inaccurate solutions. There- 
fore, we had to invert the amplitude and phase self-calibration 
after subtracting this source from the field. We used the casa 
task invgain (Hales 2016) to invert the calibration table and 
applied it to the self-calibrated visibilities. The final image ob- 
tained had a beamwidth size of 3.6" x2.1", and an rms noise 
level of 24 uJy beam !. 

To extract the light curve, we imaged the field using a 60- 
s time interval and a mask around the target with the cAsA 
tclean task, setting niters-10 and cycleniter=10. We se- 
lected these values after evaluating random images where the 
target was either detected or not detected, and ensuring that this 
choice did not result in false positive detections. Then, we used 
the imfit task to obtain the peak flux density of the target and 
the check source, J102348.2+004017, for each 60-s image to in- 
vestigate possible systematic effects. To measure the rms noise 
level around J1023 and check source, we placed four boxes 
around each source and used the imstat task to average the 
rms noise level within each box. To build the light curve, we 
considered a source as detected whenever the flux density value 
obtained with imfit was >30; otherwise, we took 3 times the 
rms noise level as an upper limit. To take into account possi- 
ble systematics affecting the field of view, we investigated pos- 
sible correlations between the light curves of the target and the 
check source. No evidence of such correlations was found. Fi- 
nally, we further flagged bad 60-s images and removed intervals 
with odd results, such as sudden increases in the flux density of 
J102348.2+004017. 


2.6. FAST (Night 1) 


FAST (Jiang et al. 2019, 2020; Qian et al. 2020) observed J1023 
starting on 2021 June 4 at 08:47:45 UTC. The total integration 
time was 3 hr (the first two minutes and the last two minutes were 
dedicated to the injection of an electronic noise signal, referred 
to as a CAL). The observation covered a range of orbital phases 
from approximately 0.5 to 1.1, thus encompassing the epoch of 
passage of the pulsar at the inferior conjunction of the orbit. The 
observation sampled a total of 21 low mode episodes according 
to the simultaneous XMM-Newton observation (see Fig. 1). The 
central frequency was 1.25 GHz, with a range from 1.05 GHz 
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to 1.45 GHz, including a 20 MHz band edge on each side. The 
average system temperature was 25 K. The recorded FAST data 
stream for pulsar observations is a time series of total power per 
frequency channel, stored in PsnriTS format (Hotan et al. 2004) 
from a ROACH-2-based backend^, which produces 8-bit sam- 
pled data over 4k frequency channels at a cadence of 49.152 us. 

We searched for pulsations with either a dispersion signature 
or instrumental saturation in all data collected during the obser- 
vation. We performed four different types of data processing: 

(1) Data folding based on the ephemeris derived from the X- 
ray data: we analysed the data using the ephemeris derived from 
the strictly simultaneous X-ray data collected by XMM-Newton 
(Miraval Zanon et al. 2022). We separately folded the FAST data 
from the entire observation and those obtained by stacking the 
data acquired during all the periods of high X-ray mode and all 
the periods of low X-ray mode. However, no radio pulses were 
detected with SNR>6 in either case. 

(ii) Dedicated search for persistent periodic radio pulses: for 
the sake of robustness, we created de-dispersed time series for 
each pseudo-pointing over a Dispersion Measure (DM) range 
between 3 pc cm^? and 100 pc cm^?. This range covers all uncer- 
tainties of the semi-blind search, as J1023 was discovered at a 
DM of 14.3 pc cm? (Archibald et al. 2009). We set the step size 
between subsequent trial DMs equal to the intra-channel smear- 
ing over the entire band. This ensures that any trial DM deviat- 
ing from the source DM by ADM will not cause extra-smearing 
beyond the intra-channel smearing. We used the PulsaR Explo- 
ration and Search TOolkit pipeline (PREsTo; Ransom 2001) to 
search for periodic signals in the power spectrum for each trial 
DM. We allowed the power of the highest harmonic component 
of the signal to drift in the Fourier domain up to a maximum 
number of 200 frequency bins and 50 frequency derivative bins 
(jerk search; Andersen & Ransom 2018; see also Wang et al. 
2021). This corresponds to an acceleration of 90.93 m s? and a 
jerk of 0.002 m s^? for the fundamental component. After check- 
ing all signal candidates with SNR»5, we identified them as 
narrow-band radio frequency interference (RFI). 

(iii) Search for single pulses: we used the dedicated search 
scheme in (ii) to de-disperse the data. We then applied 14 box-car 
width-match filter grids distributed in logarithmic space from 0.1 
to 30 ms. A zero-DM matched filter was used to mitigate RFI in 
the blind search. All the resulting candidate plots were visually 
inspected, but most were determined to be RFI. No pulsed radio 
emission with a dispersive signature was detected with SNR>5. 

(iv) Search for baseline saturation: we acknowledge that 
FAST would become saturated if the radio flux density reached 
levels as high as kJy to MJy. Therefore, we also searched for 
saturation signals in the recorded data-set. We looked for epochs 
where at least half of the channels satisfied one of the follow- 
ing conditions: 1) the channel was saturated (255 value in 8-bit 
channels), 2) the channel had a zero value, or 3) the rms devi- 
ation of the bandpass was «2. We did not find any instances of 
saturation lasting 20.5 s during the observational campaign. 

We then calibrated the noise level of the baseline with noise 
CAL injection, and derived a 5-c flux density upper limit for per- 
sistent radio pulsations of 2.0+0.5 uJy at 1.25 GHz considering 
the entire integration time of 2.8 hr, assuming a pulse duty cy- 
cle between 0.05—0.3 (which is typical for MSPs; see the ATNF 
pulsar catalogue; Manchester et al. 2005) and disregarding the 
possible scintillation effect. Within the uncertainty of the flux 
measurements and stability threshold, no flux variations were 


^ https://casper.astro.berkeley.edu/wiki/ROACH-2_ 
Revision 2 
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seen during the observation, consistent with an apparent lack of 
scintillation. The 5-c flux upper limits on the pulsed flux den- 
sity of J1023 during the high and low modes can be evaluated 
by scaling the above-mentioned limit for the specific integration 
time-spans during the high and low-mode periods in the X-rays: 
we obtain upper limits of 2.2+0.6 wy in the high mode (2.3 hr) 
and 4.7+1.0 uJy in the low mode (0.5 hr). 


2.7. NICER/XTI observations (Night 2) 


The X-ray Timing Instrument (XTI) mounted on NICER (Gen- 
dreau et al. 2012) observed J1023 between 2021 June 26 at 
21:56:04 UTC and 2021 June 27 at 01:33:51 UTC, for an ex- 
posure time of ~6.1 ks. The raw, unfiltered data were processed 
and analysed using tools within the NiceRDAs software package 
and the latest version of the calibration files at the time of the 
analysis (cALDB v. 20221001)°. We extracted the background- 
subtracted light curve in the 0.5-10 keV energy range with a 
time bin of 10s using the NIicERL3-Lc high-level pipeline, which 
estimates the background level using the “Space Weather” back- 
ground model. 


2.8. Swift/XRT observations (Nights 1 and 2) 


The Swift (Gehrels et al. 2004) X-ray Telescope (XRT; Burrows 
et al. 2005) observed J1023 in the photon-counting mode (time 
resolution of 2.5s) both during Night 1 and Night2. The first 
observation started on 2021 June 3 at 18:29:46 UTC and con- 
sisted of 11 snapshots, with a total elapsed time of 251.3 ks and 
an effective exposure time of ~17.9ks. The second observation 
started on 2021 June 26 at 22:19:35 UTC and consisted of three 
snapshots, resulting in a total elapsed time of ~7.3ks and an 
effective exposure time of ~2.4ks. The raw, level-1 data were 
processed and screened using standard criteria. Source photons 
were collected from a circular aperture of radius 47.2” centred 
on J1023 and background photons were collected from an an- 
nular region with inner and outer radii of 94.4” and 188.8’, re- 
spectively, also centred on J1023. We extracted a background- 
subtracted time series with a time bin of 100s, suitable to sample 
adequately the X-ray mode switching (see Figs. | and 3). 


2.9. Swift/UVOT observations (Nights 1 and 2) 


The Swift Ultra- Violet/Optical Telescope (UVOT; Roming et al. 
2005) observed J1023 simultaneously with the XRT. The ob- 
servations were performed in both image and event modes us- 
ing the UVM2 filter (central wavelength 2260 A; full-width at 
half-maximum 527 A). Data were processed and analysed in the 
standard way. For the source photons, we used a circular ex- 
traction region of radius of 5’’ centred on the source position. 
For the background photons, we adopted a circular extraction 
region of radius 20” away from source. For the data in image 
mode, photometry was carried out using the UvorsouRcE task. 
For the data in event mode (time resolution of 11 ms), we used 
the COORDINATOR tool to convert raw coordinates to detector and 
sky coordinates, UVOTSCREEN to filter out hot CCD pixels and ob- 
tain a cleaned event list, and vvorEvrLc to extract a background- 
subtracted light curve with a time bin of 120s, which provided a 
sufficiently high SNR to detect any changes in the UV brightness 
at least during the longest periods of low X-ray mode. 


5 We followed the standard procedures outlined on the NICER 
data analysis webpage (https://heasarc.gsfc.nasa.gov/docs/ 
nicer/analysis. threads/). 


2.10. REM photometric observations (Night 2) 


We observed J1023 during Night2 with the 60-cm REM tele- 
scope (Zerbi et al. 2001; Covino et al. 2004) located at La 
Silla Observatory (Chile). We obtained 48 strictly simultane- 
ous exposures with the ROSS2 imager, each with an integra- 
tion time of 60s, with SDSS griz optical filters (Table 1). We 
reduced the images following standard procedures (bias subtrac- 
tion and division by a normalized flat-field), and performed aper- 
ture photometry with the PHOT tool in IRAF using an aperture 
size of 6 pixels. We performed flux calibration using stars in 
the field from the American Association of Variable Star Ob- 
servers (AAVSO) Photometric All-Sky Survey (APASS) cata- 
logue? (Henden 2019). We obtained simultaneous K-band ob- 
servations with the REMIR instrument. A total of 300 15-s in- 
tegration images were acquired, which we combined 5 by 5 to 
correct for the background contribution. Aperture photometry 
was performed with IRAF as described above. The images were 
flux-calibrated against a selection of stars in the field from the 
2MASS point source catalogue". 


2.11. VLT/FORS2 polarimetric observations (Night 2) 


We observed J1023 with the FOcal Reducer/low dispersion 
Spectrograph 2 (FORS2; Appenzeller et al. 1998) mounted on 
the 8.2-m VLT at Paranal Observatory (Chile) in polarimet- 
ric mode using the optical filter RLS PECIAL + 76 (R; cen- 
tral wavelength 655 nm, full-width-half-maximum 165 nm) dur- 
ing Night 2. The observations were performed under thin cirrus 
clouds conditions, with seeing varying from 0.4" to 0.9". A total 
of 156 short (10 s) exposures were taken, resulting in a total ob- 
serving time of ~2.41 hr on source, overheads included (~50% 
of the 4.75 hr orbital period). 

A Wollaston prism was inserted in the optical light path to 
split the incident radiation into two orthogonally polarised light 
beams called ordinary (o-) and extraordinary (e-) beams. A Wol- 
laston mask was introduced to avoid overlap of the beams on 
the CCD. In addition, a rotating half-wave plate (HWP) was in- 
serted, allowing images to be taken at four different angles (®,) 
of the plate with respect to the telescope axis: ®; = 22.5°(i — 
1) i = 1,2,3,4. This step is essential to obtain a polarisation 
measurement, since the set of images at the four different an- 
gles must be combined to evaluate the level of linear polarisa- 
tion (as described below). Therefore, a total of 39 polarisation 
epochs were acquired with the program. All images were pro- 
cessed by subtracting an average bias frame and dividing by a 
normalized flat field. Aperture photometry was performed using 
the daophot tool (Stetson 19872), using an aperture of 6 pixels. 
The normalised Stokes parameters Q and U for the linear polar- 
isation (LP) were evaluated following the methods described by 
Baglio et al. (2020) and references therein. These values were 
then corrected for instrumental polarisation by observing a non- 
polarised standard star? (WD 1620-391) using the same setup 
during the same night. A polarisation level of « 0.296 was mea- 
sured, which is consistent with the typical instrumental polari- 
sation level expected for the FORS2 instrument as reported on 
the ESO website? («0.396 in all optical bands). The Stokes pa- 


6 http: //www.aavso.org/download-apass-data 

7 https://irsa.ipac.caltech.edu/applications/2MASS/IM/ 
8 https://www.eso.org/sci/facilities/paranal/ 
instruments/fors/tools/FORS1 Std.html 

3 http://www.eso.org/observing/dfo/quality/FORS2/ 
reports/FULL/trend report PMOS inst. pol FULL.html 
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rameters evaluated for the non-polarised standard star were then 
subtracted from the Stokes parameters of J1023 at each epoch. 
We then used an algorithm to estimate the level and angle 
of linear polarisation. This algorithm is based on the evalua- 
tion of the quantity S (0) for each of the HWP angles (e.g. di 
Serego Alighieri 1998 and references therein; see also Covino 
et al. 1999; Baglio et al. 2020). S(®) is the component of the 
Stokes vector that describes the LP along the ® direction (e.g., 
f*(9)/ FD) 1) (Amie 


S(0°) = Q). It is defined as: 
=i 
— + J A (1) 
u (®)/ fu (0) u (®)/ fu (D) 


where © is the angle of the HWP, f°(®) and f°(®) are the or- 
dinary and extraordinary fluxes for a given angle ®, and f? (®) 
and f*(®) are the same quantities evaluated for the non-polarised 
standard star. S(®) is related to the polarisation level P and 
angle 0 by S(®) = P cos2(0 — 4). The parameters P and 0 
were estimated by maximizing the Gaussian likelihood func- 
tion and using this as the starting point for a Markov Chain 
Monte Carlo (MCMC) procedure (Hogg & Foreman-Mackey 
2018) (see Baglio et al. 2020 for further details on the algorithm). 
The best fit parameters were then calculated as the median of the 
marginal posterior distributions of these parameters, while the 
uncertainties were estimated as the 16—84th percentiles of the 
distributions. For non-detections, the 99.97% percentile of the 
posterior distribution of the parameter P was used to calculate 
an upper limit. 

After fitting the S () curve, we obtained an LP value that is 
already corrected for instrumental effects. However, the polarisa- 
tion angle still needs to be corrected using the Stokes parameters 
of a polarised standard star. For this purpose, we observed the 
star Hiltner 652 and measured its Stokes parameters to determine 
its polarisation angle. This was then compared to the tabulated 
value, and a difference of 1.71? + 0.30? was applied to the polar- 
isation angle of J1023. However, the resulting linear polarisation 
will still include an interstellar component. Evaluating the inter- 
stellar contribution is not straightforward. In principle, we could 
use a group of reference stars in the same field as the target in 
Eq. 1, under the assumption that these stars are non-polarised. 
This way, the LP of the target would be corrected for both instru- 
mental and interstellar contributions. However, it is important to 
ensure that the chosen reference stars all have the same level of 
polarisation, which would be a good indication of their intrinsic 
non-polarised nature (see e.g. Baglio et al. 2020). Unfortunately, 
not only is the field of view of J1023 sparse in the optical, but 
the very short exposure times (10s) exclude all faint (R > 19) 
stars from this study. We are therefore left with only one possi- 
ble comparison star in the field, which makes it difficult to verify 
whether the chosen star is indeed intrinsically non-polarised. As 
a result, we used the fluxes of the non-polarised standard star in 
Eq. 1. The field star will be used as a comparison for the analysis. 

We can evaluate the maximum interstellar contribution to 
the LP of a source, referred to as Pin, through considerations 
of its absorption coefficient in the V-band, Ay: Pin(%) < 3Ay 
(Serkowski et al. 1975). For J1023, the best-fitting absorption 
column density, Ny, according to XMM-Newton observations is 
5.0*05 x 107° cm? (Campana et al. 2019). Using the relation be- 
tween Ny and Ay in our Galaxy from Foight et al. (2016), we 
calculate Ay = 0.17 +0.01 mag and therefore Pin < 0.52%. This 
is very low and comparable to the instrumental polarisation mea- 
sured for the instrument configuration. It is also comparable to 
the upper limit of the interstellar polarisation measured for J1023 
(Baglio et al. 2016). 


se»-[ 
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Table 2. ALMA continuum maps 


Weighting Synthesized beam PA rms* 
("x") C) (uly beam") 

Natural 0.87 x 0.53 —74.5 6.1 

Briggs, robust=0 0.55 x 0.30 —87.3 10.6 


* Rms noise level at the phase center. 


2.12. ALMA observations (Night 2) 


The ALMA observations were conducted using the 12 m ar- 
ray and 38 antennas in the C43-7 configuration (baselines 22— 
2492 m), with a field of view (FoV) of 59.72” and a maximum 
recoverable scale of 11.5’’. The total integration time on source 
was ~2 hr and the phase center was R.A. = 10°23™47°68; Dec. 
= +00°38’41.”00 (ICRS). The weather conditions were good 
and stable, with an average precipitable water vapor of ~1.2 mm 
and an average system temperature of ~65 K. The band-pass 
and flux calibrators used were J1058+0133 and J1256—-0547, 
while J1028—0236 was used as the phase calibrator. The con- 
tinuum image was obtained by averaging the four base-bands of 
width 1.875 GHz, each centered at 90.5003 GHz, 92.5003 GHz, 
102.5003 GHz, and 104.5003 GHz, covering a total bandwidth 
of 7.5 GHz. No spectral lines were detected within the observed 
frequency range above the 3c noise level per spectral channel of 
~0.48 mJy beam '!. 

Data reduction was carried out following standard proce- 
dures using the ALMA pipeline within casa (v. 6.1.1-15)!°. Con- 
tinuum images were produced using TCLEAN interactively and ap- 
plying a manually selected mask. These images were corrected 
for primary beam attenuation and were produced using both nat- 
ural weighting and a Briggs robust=0 weighting (Briggs 1995). 
The former maximises the sensitivity whereas the latter provides 
higher angular resolution. The mean frequency of the contin- 
uum image is 97.51 GHz, which corresponds to a wavelength 
of 3.1 mm. Table 2 lists the synthesized beam and the root mean 
square (rms) noise of each image. 


3. Results 
3.1. X-ray and UV emissions 


The X-ray light curve extracted from XMM-Newton data during 
the first night reveals flaring activity over a time span of ~1.5 hr 
at the beginning of the observation. A total of 288 switchings 
between high and low modes were observed during the observa- 
tion (see Fig. 1). Based on our thresholds for the EPIC net count 
rates in the different modes (see Miraval Zanon et al. 2022), we 
estimated that J1023 spent ~65, 17 and 1% of the time in the 
high, low and flaring modes, respectively, and the rest of the time 
switching between modes. The good time intervals (GTIs) cor- 
responding to the high and low modes were then used to extract 
the background-subtracted spectra, response matrices, and ancil- 
lary files separately for each mode. The NuSTAR light curve is 
characterised by a lower SNR than the XMM-Newton/EPIC light 
curve. Nevertheless, it is possible to recognise the same flaring 
episodes and mode switchings detected using XMM-Newton (see 
Fig. 1). We used the GTIs for the high and low modes selected by 
the strictly simultaneous XMM-Newton/EPIC light curves to ex- 
tract the background-subtracted spectra, response matrices, and 
ancillary files separately for each mode. 

The X-ray light light curve extracted from NICER data dur- 
ing the second night, shown in the top panel of Fig.3, re- 


10 https://casa.nrao. edu 
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veals 22 instances of X-ray mode switching, without any ob- 
vious flaring episodes. We applied intensity thresholds to the 
light curve to select the time intervals corresponding to the high 
(22.5 counts s^!) and low («2.5 counts s^!) X-ray modes. Based 
on this selection, we estimate that J1023 spent ~70% of the time 
in the high mode and ~30% of the time in the low mode. We also 
detected two periods of X-ray low mode in the Swift/XRT data, 
with the former also covered by the NICER observations (see the 
second panel of Fig. 3). 

The UV light curves exhibit substantial variability. Specif- 
ically, the HST light curve during Night 1 displays flaring ac- 
tivity as well as dips in the count rate that coincide with the 
low-mode episodes observed in X-rays (Fig. 1; see also Miraval 
Zanon et al. 2022; Jaodand et al. 2021). A dip in the UV bright- 
ness by ~0.7 mag is also clearly visible in the Swift/UVOT data 
during the third snapshot of Night 2. The decrease occurs over a 
span of several minutes during an episode of low mode (Fig. 3). 


3.2. Optical/NIR photometric properties 


During the first night, the source is brighter at the beginning than 
in the rest of the observation (reaching B ~ 17 — 16.5 mag), dis- 
playing a smooth quasi-sinusoidal modulation at the orbital pe- 
riod with a semi-amplitude of ~0.2 mag around an average value 
of B ~ 17.3 mag. We fit a sinusoidal function to the light curve 
by fixing the period and reference phase of the sinusoid to the bi- 
nary orbital period and the phase value obtained from the time of 
passage of the pulsar through the ascending node, respectively, 
as derived by Miraval Zanon et al. (2022). Then, we subtracted 
the best-fitting function from the observed count rates. The re- 
sulting detrended light curve is shown in Fig. 1. After initial flar- 
ing activity resembling the X-ray behaviour, the optical intensity 
appears to be flickering, but without any sign of mode switching. 
Fig. 1 also shows the NIR light curve. Variability in the form 
of dips and flickering is evident. The fractional rms magnitude 
deviation of the light curve, computed following the method 
described by Vaughan et al. (2003), is (10.1+0.2)%. No NIR 
flare was observed at the switch from the low to the high X-ray 
mode, as previously reported (Papitto et al. 2019; Baglio et al. 
2019). Furthermore, we did not detect any signs of mode switch- 
ing in the light curve. The average magnitudes measured sepa- 
rately during the high and low modes are compatible within 1c 
(J = 15.39 + 0.16 mag and J = 15.37 + 0.13 mag in the high and 
low modes, respectively), indicating that there is no significant 
correlation between the NIR variability and mode switching. 
During Night2, the optical light curves show clear variabil- 
ity on short timescales (Fig.3, fourth panel). Unlike what was 
seen during Night 1, no modulation at the system orbital period 
is detected. The average magnitudes are g = 16.65 + 0.02 mag, 
r = 16.45x0.02 mag, i = 16.20+0.03 mag, z = 16.09 0.03 mag 
(in the AB photometric system), all marginally consistent with 
those reported by Coti Zelati et al. (2014); Baglio et al. (2016). 
To estimate the amount of intrinsic variability, we evaluated the 
fractional rms magnitude deviation of the light curves follow- 
ing Vaughan et al. (1994). We find a low-significance fractional 
rms of 6.2 + 3.4%, 6.9 + 2.3% and 12.7 + 4.7% at a time res- 
olution of 160s in the g, r and z bands, respectively, while no 
excess intrinsic variability could be detected in the i band. This 
fractional rms is comparable to what has been observed for a 
few black-hole X-ray binaries in the optical band on a similar 
timescale, such as GX 339-4 (Gandhi 2009; Gandhi et al. 2010) 
and MAXI J1535—571 (Baglio et al. 2018). In these systems, the 
variability has been attributed to the emission of a flickering jet 
that contributes up to optical frequencies. In contrast, when ac- 


cretion dominates, lower values of fractional rms are typically 
measured, as seen in the case of the black-hole X-ray binary 
GRS 1716-249 (Saikia et al. 2022). 

Unfortunately, most of the K-band images are of bad qual- 
ity, resulting in an empty field of view. We stacked together 
three of the combined images in which the field was visible 
and J1023 was detected (so a total of 15x15s integration im- 
ages, all taken during the high mode), and performed the cali- 
bration against a selection of stars in the field from the 2MASS 
point source catalogue. We obtained an average magnitude of 
K = 14.41 + 0.29 mag, which is consistent with previous values 
reported in the literature (see e.g. Coti Zelati et al. 2014; Baglio 
et al. 2016). 


3.3. Optical polarimetric properties 


Fig. 3 shows the evolution of the optical LP level during the sec- 
ond night of observations. To improve the SNR of the obser- 
vations without losing too much information on the short-term 
variability of the LP, we combined the measurements so that 
each LP data point covers 80s (4 HWP angles x 10s of inte- 
gration x 2 sets). We define an LP measurement as a detection 
when it reaches a confidence level of greater than 20. Using this 
threshold, we obtain a total of 8 LP detections and 11 upper 
limits at a 99.97% confidence level (see Table 3). The average 
level of linear polarisation, calculated as the average of all 20- 
significant detections, is 0.92 + 0.39%. This is consistent with 
previous results in the same band (R-band; Baglio et al. 2016). 
The polarisation angle is well constrained at a value of (6+18)° 
(corrected for the polarisation angle of the standard star). 

When detected, the LP appears fairly constant throughout the 
observations. We investigated any potential changes in the LP 
level as a function of the X-ray mode (low/high) by averaging 
all images taken during the high and low modes. We found that 
the LP was 0.65 + 0.43% (with a 99.97% confidence level upper 
limit of 2.12%) in the high mode and 0.36+9:31% (with a 99.97% 
confidence level upper limit of 1.58%) in the low mode. 

The measured levels of LP are not corrected for interstel- 
lar contribution, which is expected to be low according to 
Serkowski's law (<0.52%; Serkowski et al. 1975). Fig. 5 shows 
the Q and U Stokes parameters on the Q — U plane for the whole 
dataset, including J1023 (black dots), a comparison star in the 
field (red squares), and the non-polarised standard star (green 
star). The non-polarised standard star is located very close to 
the origin of the Q — U plane, indicating that the instrumen- 
tal polarisation of VLT/FORS2 is low, as expected. The Stokes 
parameters of the comparison star, on the other hand, are scat- 
tered around a common point, consistent with the position of 
the non-polarised standard within ~1.56, where 6 is the standard 
deviation of the points. Any deviation from the position of the 
non-polarised standard on the plane is a good estimate of the 
interstellar polarisation of the field (assuming that the compari- 
son star is intrinsically non-polarised). Interestingly, the Stokes 
parameters of J1023 are clustered around a common point in a 
different region of the Q — U plane, with a higher scatter than the 
comparison star due to larger uncertainties. For the two time- 
series of the Q and U Stokes parameters, we are unable to claim 
any detection of short-timescale variability. 


3.4. Millimetre emission 


Fig.6 shows the ALMA 3.1 mm continuum emission towards 
J1023. Two sources were detected within the FoV of ALMA. 
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Fig. 1. Temporal evolution of the X-ray, UV, optical, NIR and radio emissions of J1023 on 2021 June 3—4. The light curves are shown in decreasing 
order of energy band from top to bottom with the XMM-Newton/EPIC one as the reference for all others (NuSTAR, HST/STIS, XMM-Newton/OM, 


NTT/SOFL VLA). The XMM-Newton X-ray light curve is shown over 


the time interval in which the three EPIC instruments collected data 


simultaneously. The OM residual count rates were obtained after removing the orbital modulation from the time series (see the text for details). 
Only VLA detections are shown for displaying purpose. The dark yellow segment in the bottom panel indicates the time range covered by the FAST 
observation. The orange shaded areas mark the time intervals of the low X-ray mode detected by XMM-Newton/EPIC covered by observations 


with other instruments. 


J1023 was detected at a significance level of ~15c. The other 
source, which we designate ALMA J102349.1+003824 follow- 
ing the naming convention of new sources discovered by ALMA, 
is located ~22” south-east of J1023 and was detected at a signif- 
icance level of —^9c. The panels on the right of the figure dis- 
play the zoom-in maps around each source, with the continuum 
image obtained with natural weighting displayed in color scale 
and the black contours derived from the image with robust-8. 
Table 4 reports the peak coordinates, the peak intensity, the in- 
tegrated flux, and the deconvolved size obtained by fitting a 
two-dimensional Gaussian function to the high-resolution image 
(robust-6). 


The 3.1-mm continuum emission from J1023 has a compact 
shape with some elongation towards the southeast (shown in the 
top right panel of Fig. 6 using black contours). The millimetre 
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continuum emission peaks ~0.32” south of the optical position. 
The deconvolved size obtained from the 2D Gaussian fit in the 
robust=0 image is 0.305" x 0.205”, which corresponds to ap- 
proximately 418AU x 281AU (~6.3x 10? cm x 4.2x10!> cm) at 
a distance of 1.37 kpc (this is approximately 7x10* times the or- 
bital separation of the system; Archibald et al. 2013). The emis- 
sion from ALMA J102349.1 003824 also appears compact (see 
the bottom right panel of Fig. 6). The 2D Gaussian fit yields a 
point-like source with a FWHM as large as 0.39" x0.22" (see Ta- 
ble 4). The results of the searches for counterparts of this source 
at other wavelengths are reported in Appendix A. 


To determine the spectral index of J1023 within the ALMA 
bandwidth, we imaged the continuum emission from the 
lower and upper sidebands, which have central frequencies of 
91.47886 GHz and 103.5115 GHz, respectively. The resulting 
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Fig. 2. Zoom-in of the time series collected on Night 1 over the time span of maximum overlap among the different telescopes. Arrows indicate 


3-c upper limits. 


spectral index is —0.1 +0.8, which is consistent with previous re- 
ports for the average cm emission (Deller et al. 2015; Bogdanov 
et al. 2018). 


To probe variability in the millimetre emission, we obtained 
images scan by scan (typically with a duration of ~4—5 minutes), 
using a natural weighting. We then computed the peak intensi- 
ties and flux densities considering the emission above 30 level. 
The bottom panel of Fig.3 shows the evolution of the flux den- 
sity over time. This quantity was found to vary by a factor of ~5 
within the range ~30-160 uJy on a timescale of about 5 minutes. 
In some cases, the flux density remained below the detection sen- 
sitivity, resulting in upper limits of the order of a few tens wy. 
Fig. 3 also shows the ALMA time series extracted using time in- 
tervals of variable length corresponding to the durations of the 
X-ray-mode episodes. Two flaring episodes can be clearly seen. 
We measured average flux densities of 89x13 Jy in the high 
mode and 134+23 uJy in the low mode, excluding the two flar- 
ing episodes (see Table 5). Overall, the anticorrelated variability 


pattern appears to be more pronounced between the cm and X- 
ray emission than between the millimetre and X-ray emissions. 


3.4.1. Millimetre flares 


The integration times adopted above are effective for studying 
the variability of the millimetre emission on minute timescales, 
but may average out large enhancements of the millimetre flux 
on much shorter timescales (of the order of seconds). Therefore, 
we performed a dedicated search for short flaring episodes by 
extracting images of different durations using the entire dataset. 
This search resulted in the detection of the two above-mentioned 
episodes of bright millimetre emission (Fig. 4). No other remark- 
able events were detected. The first episode (hereafter "Flare 
1”) occurred on 2021 June 26 within the time range 22:19:29 
— 22:19:44 (UTC) at the switch from high mode to low mode 
and reached a peak flux density of 1.01+0.08 mJy. The second 
event (“Flare 2”) took place during the time range 23:48:29 — 
23:48:44 (UTC) and reached a peak flux density of 0.7+0.1 mJy. 
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Fig. 3. Temporal evolution of the X-ray, UV, optical and millimetre emissions and of the optical linear polarisation of J1023 on 2021 June 26- 
27. The light curves are shown in decreasing order of energy band from top to bottom (NICER, Swift/XRT, Swifi/UVOT, REM, VLT/FORS2 
polarisation, ALMA). In the bottom panel, the ALMA time series are plotted using either bins of length 4-5 min (cyan) or bins of variable length 
corresponding to the duration of the X-ray mode episodes identified by the strictly simultaneous NICER and Swift data (blue). In some panels, 
the sizes of the markers are bigger than the uncertainties and/or the time bins. Arrows indicate 30 upper limits. The orange shaded areas mark 
the time intervals of low X-ray mode detected by NICER and Swift/XRT covered by observations with other instruments. A plot of the X-ray and 
millimetre light curves around the epoch of the millimetre flares is shown Fig. 4. 


The lack of continuous coverage in the millimetre band before 
and during the switch from the high to the low mode does not 
allow us to determine the exact epoch of the flare onset. We note 
that these values refer to the flux at the peak of the flares and are 
therefore higher than those obtained in the X-ray mode-resolved 
analysis (see Fig. 3 and Table 6), where the emission is averaged 
over longer time intervals. We observed a significant increase 
of the peak intensity above the average level: a factor of ~4 for 
Flare 1 and a factor of ~7 for Flare 2. 


We observed a clear change in the morphology of the mil- 
limetre continuum emission during the two flaring episodes (see 
white contours in Fig. 7). The two-dimensional Gaussian fitting 
(see Table 6) indicates that the emission is elongated along the 
northeast to southwest direction during Flare 1 (roughly perpen- 
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dicular to the ALMA 3.1 mm emission obtained with all the vis- 
ibilities), while it is marginally resolved in only one direction 
during Flare 2. The positional accuracy of ALMA images can be 
estimated as beampwum X (0.9 x SNR)! !!, where beamgwum 
is the full-width-half-maximum synthesized beam in arcsec and 
SNR is the signal-to-noise ratio of the image. The positional ac- 
curacy for Flare 1 is ~ 0.12—0.15", while for Flare 2 it is slightly 
smaller at ~ 0.1". However, the actual absolute positional ac- 
curacy may be worse than the nominal values, potentially by a 
factor of two or more, depending on the atmospheric phase con- 
ditions during the observation. The offset between the ALMA 
peak positions and the radio position (Deller et al. 2012) is 20.4" 


!! https://almascience.eso. org/documents- and-tools/ 
cycle9/alma-technical-handbook 
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Fig. 4. X-ray (top) and millimetre (bottom) time series collected around 
the epochs when the two millimetre flares were detected. The X-ray 
light curve is binned at a time resolution of 10s while the millimetre 
light curve is binned at 15 s. Arrows indicate 3-0 upper limits. Note the 
lack of coverage in the millimetre band for Flare 2 before and during 
the switch from the high to the low mode in the time interval between 
two consecutive ALMA scans. Note also the different scales on the hor- 
izontal axis in the panels on the left and on the right. 
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Fig. 5. Stokes Q and U measurements for the whole VLT/FORS2 
dataset. Black dots, red squares and the green star indicate the Stokes 
parameters of J1023, a comparison star in the field and the non-polarised 
standard star, respectively. 


for Flare 1 and ~0.3” for Flare 2. Given the stable weather con- 
ditions during the observation, it is possible that this offset is 
genuine for both flares. However, caution must be exercised, as 
this offset is comparable to the positional accuracy of ALMA. 


Table 3. Linear polarisation detections and angles with a significance 
>2ø. For non-detections, upper limits are quoted at a confidence level 
of 99.97%. The first column indicates, in parentheses, which mode the 
MJD refers to (‘L’ for low mode, ‘H’ for high mode, ‘?’ when unknown 
due to the lack of simultaneous X-ray observations). 


Epoch (MJD) LP (95) PA() Upper limit (%) 
59391.95402 (7) — 07605 273250 
59301.05017(7) 0.814838 — -3.5055 

59391.96436 (?) - - 1.77 
59391.96949 () | 1.0007 . -2.92:025 

59391.97459 (?) —— 1.04*022 1569 

59391.97971 (?) - - 1.99 
59391.98485 (L) - - 1.82 
59391.99000 (H) - - 3.42 
59391.99517 (H) - - 4.77 
59392.00034 (H) - - 2.33 
59392.00549 (H) 1.12+0.32 21.96+814 

59392.01063 (?) - - 2.02 
59392.01578 (?) - - 1.22 
59392.02094 (?) — 0.63:030  —6.46+!238 
59392.02608 (?) 0.96+0.29 138723 

59392.03126 (?) - - 1.61 
59392.03643 (?) - - 1.27 
59392.04160 (?) 1.037038 9.3941) 39 

59392.04680 (H) - - 2.25 


3.5. Radio emission 
3.5.1. Continuum emission 


The radio light curve extracted from VLA data is shown in Fig. 1. 
J1023 is not detected most of the time, especially during the 
high mode intervals. However, it is detected during almost all 
episodes of low mode (here the uncertainties on the flux densi- 
ties are reported at the 1o confidence level). The average radio 
flux densities measured separately in the high and low mode are 
reported in Table 7. These values are consistent with those re- 
ported by Bogdanov et al. (2018) within the uncertainties. This 
indicates that the radio flux densities measured separately in the 
two modes are remarkably stable over time scales of years. 


3.5.2. Non-detection of pulsed emission 


The limit on the radio pulsed flux density derived from the FAST 
observations is a factor of 25000 smaller than the flux density 
measured at similar frequencies when J1023 was active as a ra- 
dio pulsar (Archibald et al. 2009) and a factor of a few tens 
more constraining than the limits previously evaluated (again at 
similar frequencies) as soon as J1023 transitioned to the sub- 
luminous X-ray state (Stappers et al. 2014; Patruno et al. 2014). 
The limit on the radio pulsed flux density in the high mode is also 
at least four orders of magnitude smaller than the value predicted 
by extrapolating the best-fitting power-law model for the SEDs 
of the optical-to-X-ray pulsed emission (Papitto et al. 2019; Mi- 
raval Zanon et al. 2022). 
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Table 4. Parameters of the sources detected at 3.1 mm with ALMA. 


Source R.A (ICRS) Dec. (ICRS) L^ Ss, Deconvolved Size* PA 

(hh:mm:ss) (Qi) (uJy beam*!) (uJy) ("x") (°) 
PSR J1023 10:23:47.689+0.001  00:38:40.677+0.004 103.7+3.6 144.6+8.0 (0.305+0.060)x(0.205+0.031) 86229 
ALMA J102349.1 . 10:23:49.131+0.002 00:38:23.72+0.02 9515 1094-29 <0.39x0.224 


a Peak intensity. " Integrated flux. © FWHM obtained from 2D Gaussian 


fits. ^ Component is a point source. 
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Fig. 6. Left: ALMA 3.1 mm continuum image showing the two detected sources, J1023 and ALMA J102349.1+003824. The white dotted boxes 
display the field of view shown in the right panels. The ALMA synthesized beam (0.87" x 0.53", P.A.=—74.5°) is shown at the bottom right 
corner. Right: Close-up images of J1023 (top) and ALMA J102349.1+003824 (bottom). In both panels, the color scale corresponds to the image 
with natural weighting with the 3c contour indicated in white, where œ is listed in Table 2. The black contours are derived from the ALMA 
image obtained with robust-0. Contours range from 3c to 9c in steps of 1c, where o = 9.1 uJy beam™! for J1023 and o = 15 uJy beam! for 
ALMA J102349.1+003824. The white cross depicts the radio position of J1023 (Deller et al. 2012). A faint emission tail is detected at 3a level 


extending south-west of J1023. The synthesized beams for the robust=0 and natural images are displayed in the bottom left and bottom right 


corners, respectively (see Table 2). In all panels, the color scale is in units of Jy beam '. 


1 


Table 5. X-ray and millimetre properties of J1023 during the periods of high and low X-ray modes. 


X-ray mode T’ Fx? L, Dy Deconvolved Size PA 
(10-? erg cm? s?!) (uJy beam?!) (uy) ("x") (°) 
High 1.66 + 0.03 13.8 + 0.2 74+6 89x13  (0.338+0.202)x(0.250+0.044) 105+8 
Low 1.70.1 2.1 + 0.2 11711 134+23 (0.384+0.245)x(0.260+0.125) 141+45 


Notes. * Best-fitting photon index 


for the absorbed power-law model. 


> Unabsorbed flux over the 0.3-10 keV energy band. 


3.6. The low-mode ingress and egress timescales 


We first examined the low mode ingress and egress times in the 
XMM-Newton observation. For each low-mode episode, we used 
a simple model for the count rate light curve consisting of four 
components: a constant level for the high mode (allowing for 
a different count rate at the entrance and at the exit of the low- 
mode episode), a linear ingress in the low mode, a constant count 
rate in the low mode, and a linear egress. This model is similar in 
some aspects to what is used to model eclipses in binary stars or 
planetary transits. While we do not expect to accurately describe 
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each low-mode episode individually, we consider the results in a 
statistical sense. 


We applied this model to each low-mode episode in the 
XMM-Newton/EPIC count rate light curve, which was binned 
at 10s, excluding the initial flaring mode of the observation. We 
performed the same procedure on the NICER data from the sec- 
ond night, using again a 10-s time bin. We fit the model to 56 
low-mode episodes for XMM-Newton and 10 for NICER. 


As shown in Fig. 8, the ingress time is shorter than the egress 
time from the low mode (~ 80%). Interestingly, we observe that 
the two low-mode instances detected in the NICER data and 
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Fig. 7. Contour levels for the millimetre emission detected during 15-s time intervals at the peak of Flare 1 (left panel) and Flare 2 (right panel) are 
shown using white solid lines overlaid on the naturally weighted ALMA 3.1 mm image extracted using the whole integration time (color scale). 
These contour levels start at 3o and increase in steps of 1c", where c is the rms of the image, ~90 Jy beam"! . Both panels display the same color 
scale range in units of Jy beam". Time intervals (in UTC) are displayed at the top of each panel. The red cross depicts the radio position of J1023 
(Deller et al. 2012). The synthesized beams, with extension 0.696” x 0.540" and position angle P.A.— —79.92° for Flare 1 and 0.931" x 0.525" 
and P.A.= —74.21° for Flare 2, are shown at the bottom-right corner of each panel. 
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Table 6. Millimetre properties of J1023 during the two flaring episodes. 


Flare R.A (ICRS) Dec. (ICRS) PES S,? Deconvolved Size PA 
(hh:mm:ss) Ce) (uJy beam!) (uJy) ("x") 0) 

Flare1 | 10:23:47.70520.006 — 00:38:40.920.11 426326 1006-84 (0.96+0.89)x(0.40+0.10) 13.3+6.9 

Flare2 — 10:23:47.702+0.002  00:38:40.63«0.01 7701262 675x114 b b 


Notes. * These values are measured at the peak of the flares over time intervals of length 15 s and are therefore higher than those obtained in the 
X-ray mode-resolved analysis and plotted in Fig. 3, where the emission is averaged over longer time intervals. ^ Source marginally resolved in 
only one direction. 


coinciding with short-duration flares detected by ALMA both 
have comparable ingress and egress times (see the left panel of 
Fig. 8). This is uncommon and may be related to the detectability 
of flares at millimetre wavelengths. The distributions of XMM- 
Newton low-mode ingress and egress times are shown in the 
right panel of Fig.8. The average logarithmic ingress time is 
(10.4+6.9) s, which is consistent with the time binning. The aver- 
age logarithmic egress time is (30.6-- 8.2) s. We also investigated 
whether there was any relationship between ingress time, egress 
time, duration of the low mode, count rate at the beginning and 
end of the low mode, and count rate in the low mode. However, 
we did not find any statistically significant correlations. 

Although the model we used to describe the low-mode 
episodes is simple, it is clear that the duration of the ingress time 
in the low mode is short and consistent with the adopted time 
binning, and the egress time from the low mode is measurable 
and longer. 


4. Discussion 
4.1. The physical picture 


Our campaign, along with previous findings, allows us to draw a 
comprehensive physical picture of the high-low mode switches 


in J1023, which involves a permanently-active rotation-powered 
pulsar, an accretion disc and episodes of discrete mass ejection 
on top of a compact jet. In the following, we outline the physi- 
cal scenario while in Sect. 4.2 we will describe in detail the SED 
modeling in support of our scenario. Fig. 9 shows a visual repre- 
sentation of our scenario. We do not discuss the flaring activity 
here as it was not covered over a wide energy range. 

We assume that during the high mode the accretion flow is 
kept just beyond the light cylinder radius (280 km) by the ra- 
diation pressure of the particle wind from an active rotation- 
powered pulsar (Papitto et al. 2019). This causes the innermost 
region of the truncated, geometrically thin accretion disc to be 
replaced by a radiatively inefficient, geometrically thick flow, as 
expected for low-accretion rates (Narayan & Yi 1994). In this 
configuration, synchrotron radiation at a shock front that forms 
where the pulsar wind interacts with the inner flow produces 
the bulk of the X-ray emission as well as X-ray, UV, and op- 
tical pulsed emission (Papitto et al. 2019; Veledina et al. 2019). 
The collision of high-energy charged particles from the pulsar 
wind with the inflowing matter in the disc increases the ionisa- 
tion level of the disc. As a result, the inflowing matter is drawn 
into the pulsar magnetic field and accelerated, producing a com- 
pact jet of plasma that streams out along the direction of the pul- 
sar magnetic field lines. According to our SED modelling, the 
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Fig. 8. Results of the statistical analysis of the X-ray mode switching timescales. Left: Duration of the low-to-high mode switch (egress) as a 
function of the duration of the high-to-low mode switch (ingress) computed using data from XMM-Newton/EPIC (in red) and NICER (green). The 
circled green points represent the values for the ingress associated with the two millimetre flares detected with ALMA and the subsequent egress. 
The black solid line indicates where the two quantities are equal. Right: distribution of the low-mode ingress (in green) and egress (in red) times 


in logarithmic scale from XMM-Newton and NICER data. 


corresponding jet synchrotron spectrum is optically thick up to 
the mid-infrared band and breaks to optically thin above a fre- 
quency of (2.441?) x 10! Hz, giving only a small contribution 
at X-ray energies (~8—18% in the 0.3-10 keV band; see Fig. 11 
and Table 8). 


The switch from the high to the low mode can be connected 
to the occurrence of short millimetre flares. There is strong evi- 
dence that at least one such flare is directly involved in triggering 
the switch (see Fig. 3). We interpret these flares as the fingerprint 
of a discrete mass ejection that marks the removal of the inner 
flow (see Sect. 4.1.1), similar to what has been observed in other 
black-hole and NS X-ray binaries (e.g. Cyg X-3, Circinus X-1, 
V404 Cygni, MAXI J1535-571, MAXI J1820+091; Egron et al. 
2021; Fender et al. 1998, 2023; Russell et al. 2019; Homan et al. 
2020), albeit in different accretion regimes. As soon as the inner 
flow is expelled in the discrete ejection, synchrotron emission 
from the shock no longer takes place, resulting in a drop of the 
X-ray flux and the quenching of X-ray, UV and optical pulsa- 
tions. The system enters the low mode. 


During the low mode, the pulsar wind can still penetrate 
the accretion disc, collide with the inflowing matter, and launch 
the compact jet, contributing to the multiband emission in the 
same way as in the high mode. The X-ray emission is likely 
to be (almost) entirely from optically thin synchrotron radiation 
from this jet (Fig. 11, right panel), with a power-law index of 
a = —0.78 + 0.04 (where the flux density F, scales with the ob- 
serving frequency v as F, œ v”; see Table 8). At the same time, 
the emission associated with the discrete ejection progressively 
shifts towards lower frequencies as the ejected material expands 
adiabatically during the low mode and quickly evolves from op- 
tically thick to optically thin (Bogdanov et al. 2018). This emis- 
sion component has a steep spectrum (a = —0.8*03) and pro- 
vides an additional flux contribution above that of the underly- 
ing optically-thick emission from the compact jet at frequencies 
below the millimetre band (see Figs. 1 and 3). 

Eventually, the flow from the disc starts to refill the inner re- 
gions just outside the NS light cylinder on a thermal timescale 
(tens of seconds) due to its advection-dominated nature. The re- 
filling of the disc re-establishes the region of shock with the pul- 
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sar wind, leading to increased emission and pulsations at X-ray, 
UV, and optical frequencies (Papitto et al. 2019; Miraval Zanon 
et al. 2022; Illiano et al. 2023). The system thus makes a switch 
back to the high mode. 


As shown in Sect. 3.6, the low mode ingress and egress times 
show a statistically significant difference, with the average egress 
time being longer than the median ingress time. In our picture, 
the ingress in the low mode is caused by the ejection of the inner- 
most regions of the disc. This will occur on a dynamic timescale. 
We can estimate the distance at which the innermost edge of the 
disc in J1023 must recede away from the NS at the high-to-low 
mode switch by assuming that most of the X-ray emission in the 
low mode is due to the jet and hence imposing that the lumi- 
nosity of the synchrotron emission at the shock front is at most, 
say, 10% of that observed in the low mode. This gives a radial 
distance of row = [Lsa/ (0.1Lx tow) ]°° ric & 20rpc «1600km 
(Lsa = 4.43 x 10*4 erg s^! is the spin-down luminosity; Archibald 
et al. 2013). This value should be considered as a lower limit: 
at larger radii, the X-ray luminosity due to the shock would be 


even smaller. The dynamic timescale is ¥r3/(GM) ~ 0.1 s at 
20 rtc. Observations from XMM-Newton and NICER show that 
current X-ray instruments are capable of probing mode switch- 
ing in J1023 on timescales as short as ~10s. In the future, pro- 
posed X-ray missions that offer instruments with much larger 
collecting areas, high throughput, and excellent time resolution, 
such as the New Advanced Telescope for High-ENergy Astro- 
physics (NewAthena; Nandra et al. 2013), the enhanced X-ray 
Timing and Polarimetry mission (eXTP; Zhang et al. 2016) and 
the Spectroscopic Time-Resolving Observatory for Broadband 
Energy X-rays (STROBE-X; Ray et al. 2019), would be able to 
sample the switch on a time scale «10 s and test our predictions. 

The ejection and replenishment of the inner flow in J1023 
during mode switches from high to low mode and from low to 
high mode are reminiscent of the ejection of the X-ray corona 
into a jet and subsequent replenishment observed in the Galactic 
microquasar GRS 1915+105 (Méndez et al. 2022). A previous 
study has demonstrated that, in the case of GRS 1915-105, the 
outer disc steadily refills the central part of the system on a vis- 
cous timescale (Belloni et al. 1997). The rise time is the duration 
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Fig. 9. Schematic visual representation of the evolution of the inner 
flow into an outflow at the switch from the high to the low X-ray emis- 
sion mode. When the system is in the high mode (top panel), a small- 
size inner flow is present, together with a faint steady jet that gives rise 
to the observed low-level radio and millimetre emission. As the pul- 
sar rotates, the pulsar wind wobbles around the equatorial plane (see 
e.g. Bogovalov 1999) and shocks off the electrons in the inner flow at 
two opposite sides (red spots) at a distance that is slightly larger than 
the light cylinder radius (~80 km). At each pulsar rotation, synchrotron 
emission at the shock at X-ray, UV and optical frequencies is modu- 
lated at the spin period at one side (bright red spot), while it is absorbed 
by material in the inner flow at the other side (light red spot) (Papitto 
et al. 2019). When the system makes a switch to the low mode (bot- 
tom panel), a bright discrete ejection is launched on top of the compact 
jet, the inner flow disappears and the shock emission is quenched. For 
displaying purposes, only the inner regions of the compact jet are repre- 
sented in both panels. The direction of the NS magnetic axis is assumed 
to be misaligned with respect to the plane of the disc. Instead, the jet is 
aligned with the spin axis, as represented in the figure. 


it takes for the heating front to propagate through the central disc. 
By assuming a subsonic velocity for the heating front, with the 
sound speed cs calculated in the cool (x 6 x 10° K) disc state, 
we can estimate the rise time to be row/cs « 20s for J1023. 


In the above picture, we have assumed that the emission from 
the compact jet remains relatively constant over time, but it is 
also possible that the jet could be partially disrupted by the dis- 
crete ejection of plasma at the switch from high to low mode. 
However, if the ejection is not energetic enough to overcome the 
magnetic fields and pressures within the jet, or if it is not dense 
enough or well-aligned to effectively interact with the jet, then 
the jet would quickly re-establish itself and appear as if it had 
never turned off. In either cases, the jet would still contribute the 
same amount of flux in both modes, especially in the radio band 
where the emission is produced at greater distances from the NS 
(on the order of light-minutes to light-hours; Chaty et al. 2011) 
and therefore any changes are on much longer timescales than 
the mode switching. 


4.1.1. Millimetre flares as the fingerprint of discrete ejecta 


Starting from the time scales of the first millimetre flare that oc- 
curred during the second night (which has more accurate start 
and end times than the second flare) and the typical time scales of 
the radio flares that occurred during the first night, we attempted 
to reproduce the light curves of similar flares by assuming they 
are caused by the motion and expansion of a blob of plasma ac- 
celerated near the pulsar. Following the work of Fender & Bright 
(2019), we assumed that a typical discrete ejection expansion is 
well-described by the van der Laan model (van der Laan 1966)". 
Since no flare in this study was observed simultaneously at two 
different frequencies, we are unable to properly constrain the 
model parameters. However, we are able to check whether we 
could reproduce flares with peak flux densities and durations that 
are consistent with what we observed. We assumed that the blob 
is launched a few hundred kilometres away from the NS with 
a mildly relativistic velocity (v ~ 0.05 c), and expands linearly 
over time. We also assumed a magnetic field that scales with the 
distance from the compact object as B œ R^!, with a magnetic 
field at the launching site of the order of 10? G (see Eq.6 by 
Papitto et al. 2019 for an expression of the post-shock magnetic 
field of an isotropic pulsar wind). Additionally, we considered an 
electron number density of 10!! 10? cm^?, which is consistent 
with the lower limit estimated from high-resolution X-ray spec- 
troscopy (Coti Zelati et al. 2018). Our calculations show that we 
can reproduce flares that last for a few seconds and reach a peak 
flux density of the order of 1 mJy at the ALMA centroid fre- 
quency. Furthermore, these flares would result in flares that last 
for a few tens of seconds and reach a peak flux density of the or- 
der of 0.2 mJy at the VLA centroid frequency. These results are 
consistent with the predictions of the van der Laan model for an 
expanding blob of plasma. 

Although we cannot firmly prove the true origin of the flares, 
we can ascribe the non-detection of flares in five out of the seven 
high-to-low-mode switches covered by ALMA to the variation 
in the intrinsic properties of the ejected material. These include 
the size, density, temperature, electron energy distribution, and 
opening angle of the ejected material. In general, smaller, less 
dense, and cooler ejecta that travel at slower bulk speeds tend 
to produce fainter millimetre emission (see e.g. Tetarenko et al. 
2017; Wood et al. 2021). Additionally, the distribution of the 
ejected material along the line of sight can also affect the de- 
tectability of millimetre flares through absorption and scattering 
effects (this is likely the main reason for the non-detection of 
pulsed radio emission as well; see Sect. 4.1.2). Overall, a combi- 
nation of these factors can explain the varying emission proper- 
ties observed at millimetre and radio frequencies during (and fol- 
lowing) different high-to-low-mode switches. These variations 
are likely to affect the time it takes for the inner flow to restore 
itself and produce the shock emission, resulting in the observed 
differences in duration among the low-mode episodes. 


4.1.2. On the non-detection of radio pulsations 


The non-detection of radio pulsations could be due to ‘intrinsic’ 
phenomena within the pulsar magnetosphere (i.e. irrespective of 
the presence of external drivers such as accreting material) such 
as intermittency and/or nulling (see e.g. Weng et al. 2022 for re- 
cent similar considerations for the system LS I +61 303). How- 
ever, these phenomena were not observed during extensive radio 


12 The calculations we made are based on van der Laan (1966), and are 
described in a Jupyter-notebook available at https: //github.com/ 
robfender/ThunderBooks/blob/master/Basics.ipynb 
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timing campaigns when the system was in the radio pulsar state 
(e.g. Archibald et al. 2013; Stappers et al. 2014), so they are 
unlikely to be the cause of the non-detection. Alternatively, the 
non-detection of pulsations may be due to the presence of matter 
enshrouding the system (see also Coti Zelati et al. 2014; Stap- 
pers et al. 2014). In this scenario, we can estimate a lower limit 
on the mass transfer rate from the companion star by requiring 
that the optical depth due to free-free absorption is larger than 
unity throughout the duration of the FAST observation. The an- 
alytical expression for the free-free optical depth is given by: 


Tg.125Ggz = 2.126 x 10?x 
" m (Mws + Mc)(X + 0.5Y) (9.157 + 21n T4) 
T3? v2P2 (Mys + 0.538Mc)3 


orb 


; Q) 


where rnc is the mass transfer rate in units of Mo yr, Mys 
and Mc are the mass of the NS and the companion star, respec- 
tively, X and Y are the mass fraction of hydrogen and helium, 
respectively, T4 is the plasma temperature in units of 10* K, vg 
is the velocity of the particle wind in units of 108 ms^! and 
Pow is the orbital period in hours (Burgay et al. 2003; Iacol- 
ina et al. 2009). We assume Mys = 1.65 Mo, Mc= 0.24 Mo, 
Po»=4.75h (Archibald et al. 2009), vg=1, X=0.7 and Y=0.3", 
and vg=1. In order to have Tg1.25GHz ~ 1, and therefore not 
detect the radio signal due to free-free scattering, we estimate 
mc z 107! Mo yr! « 6 x 10^ gs-!. This value is a factor of 
210-60 larger than the mass accretion rate estimated based on 
a systematic analysis of the X-ray aperiodic variability of J1023 
in the sub-luminous X-ray state (Linares et al. 2022). This im- 
plies that a large fraction of the accreting mass does not reach the 
compact object's surface and is instead ejected from the system. 
It is worth mentioning that the predicted mass transfer rate for 
J1023, based on the transfer of mass through the loss of orbital 
angular momentum, is  107!! Mo yr^! (Verbunt 1993). 


4.2. Spectral Energy Distributions 
4.2.1. Data preparation 


We extracted the SEDs separately in the high and low mode and 
fit them using a model that reflects the physical scenario pro- 
posed in Sect. 4.1. 

The NIR, optical and UV data obtained during Night 2 were 
de-reddened assuming Ay = 0.17+0.01 mag, calculated from the 
XMM-Newton estimate of Ny using the Ny/Ay relation reported 
by Foight et al. (2016). The absorption coefficients in the other 
bands were estimated using the equations reported by Cardelli 
et al. (1989), and are listed in the last column of Table 7. 

For the low-mode SED, we considered optical data ac- 
quired between MJD 59392.01085 and 59392.01434, totaling 
three consecutive 60-s integration images in each band; the J- 
band flux is the average of the fluxes measured with NTT/SOFI 
while the source was in the low mode, according to our XMM- 
Newton monitoring; the UV (uvm2) flux was measured on MJD 
59392.01302; radio flux densities are taken from the VLA ob- 
servation on 2021 June 3-4; the millimetre flux density was de- 
rived from stacking all the ALMA images acquired during the 


13 We note that our assumption on X and Y should be regarded as 
a zeroth-order approximation since a recent study by Shahbaz et al. 
(2022) has shown that the chemical abundances of the companion star 
of J1023 are different from those found in most X-ray binaries. How- 
ever, since ric œ (X + 0.5Y)-!, the lower limit on ric only increases by 
a factor of a few if smaller values of X and Y are assumed. 
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Table 7. Central frequencies (v+), flux densities corrected for interstel- 
lar absorption and the adopted absorption coefficients (A,) used for the 
SEDs. 


Band ve (Hz) Flux Density (mJy) A, (mag) 
Low mode High mode 
Radio 1.00x 10! 0.158 20.009 0.064 + 0.006 - 
mm 9.75x10'° 0.133 + 0.023 0.089 + 0.013 - 
K 1.40 x 104 - 1.17 x 0.32 0.02 x 0.03 
J 2.39 x 104 1.14 x 0.13 1.12 + 0.16 0.05 + 0.03 
Z 328 x 10^ 1.35 x 0.28 1.34 x 0.27 0.08 + 0.01 
i 3.93 x 104 1.32 x 0.23 1.20 x 0.21 0.11 + 0.01 
r 4.81 x 10!4 1.09 x 0.11 1.03 x 0.11 0.14 + 0.01 
g 6.29 x 10^ 0.87 + 0.10 0.84 + 0.10 0.20 + 0.01 
uvm2 — 1.34x 105^ 0.10 x 0.01 0.18 x 0.02 0.52 + 0.01 


low modes. No K-band observations have been acquired during 
the low mode. The unabsorbed X-ray flux was calculated as the 
average flux during all periods of low mode. 

For the high-mode SED, we considered K-band data ac- 
quired between MJD 59392.00237 and 59392.00289, and then 
from MJD 59392.00593 to 59392.00832 (i.e. all the images in 
which J1023 is detected, 15 x 15 s); the J-band flux is the aver- 
age of all the fluxes measured with NTT/SOFI while the source 
was in the high mode, according to our XMM-Newton monitor- 
ing; for the optical griz bands, we considered one 60-s image for 
each band, at MJD 59392.00528; the UV (uvm2) flux was ac- 
quired on MJD 59392.00413; radio flux densities are taken from 
the VLA observation on 2021 June 3-4; the millimetre flux den- 
sity was derived after stacking all the ALMA images acquired 
during the high mode episodes. The unabsorbed X-ray flux was 
calculated as the average flux during all periods of high mode. 


4.2.2. Model setup 


e Low mode. To perform the fit to the SED, we used a model 
comprising a blackbody of an irradiated star and a multi-colour 
accretion disc, together with a broken power law to describe the 
optically thick and thin components of the synchrotron emission 
from the compact jet. In addition, we considered a power law to 
model the optically thin emission from the discrete ejecta at all 
frequencies. For the stellar component, we have considered an ir- 
radiated star model (see Eqs. 8 and 9 of Chakrabarty 1998) where 
the surface temperature of the star (T) is allowed to vary and 
the distance to J1023 (D), the radius of the companion star (Re), 
the binary separation (a), the irradiation luminosity (Lir) and the 
albedo of the star (7,.) are fixed to known (or reasonable) values: 
D = 1368 pc (Deller et al. 2012), Re = 0.43 Ro (Archibald et al. 
2009), 7. = 0.1, and a = [G(Mws + MoP?, /(4ny^]! 5, where 
Mws = 1.65Mg is the NS mass, Mc = 0.24Mg is the mass of the 
companion star, Porb = 4.75 hrs is the binary orbital period and 
G is the gravitational constant. Li, was fixed to 6.5 x 10?? erg s^! , 
the value estimated by Shahbaz et al. (2019) from the difference 
of the temperatures of the illuminated and dark sides of the com- 
panion star of J1023. This value corresponds to 14% of the pul- 
sar spin-down luminosity (Archibald et al. 2013), and is of the 
same order of magnitude as the X-ray and gamma-ray luminosi- 
ties in the sub-luminous X-ray state. For the multi-color accre- 
tion disc model (Eqs. 10-15 by Chakrabarty 1998), we allowed 
the inner radius of the accretion disc where the optical/UV emis- 
sion becomes relevant (rij opyuv) to vary, while we fixed the X- 
ray albedo of the disc to 0.95 (Chakrabarty 1998) and the mass 
transfer rate to 10^! ! Mọ yr !. This is the value predicted for an 
X-ray binary system with a short orbital period hosting a main- 
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sequence star and transferring mass through loss of angular mo- 
mentum (Verbunt 1993), and is consistent with the lower limit 
estimated from the non-detection of radio pulsations (Sect. 2.6). 
The broken power-law component has the normalisation Fo, the 
two slopes œ; and o» and the break frequency Vpreak as free pa- 
rameters; the power law associated with the emission from the 
discrete ejecta has the normalisation Fi, 10 (calculated at a fre- 
quency of 10 GHz) and the slope ai, as free parameters. 

e High mode. We used a more complicated model following 
the physical picture explained above. In particular, we consid- 
ered that the X-rays are produced by the sum of the optically 
thin component of the synchrotron emission from the base of 
the compact jet and the synchrotron emission from the region of 
shock between the pulsar wind and the inner flow. For the latter, 
we used the analytical expression for the synchrotron spectrum 
by Chevalier (1998) (see his Eq. 1). For the emission at lower 
frequencies, we included the blackbody of the irradiated star 
(fixing Li, to 6.5 x 1033 erg s !), the multi-colour blackbody 
of the accretion disc, and the optically thick component of the 
synchrotron emission from the compact jet. The parameters 
of the fit are Ts, voreak, @1, @2, Fo, Tin,opt/UV> Vsync» C sync and 
Fyne, Where Vsync, @syncs Fsync are the peak frequency of the 
synchrotron emission, the slope of the optically thin synchrotron 
emission and the normalisation of the synchrotron emission, 
respectively. 


Employing both data sets (high and low mode), we per- 
formed an MCMC sampling of the posterior probability den- 
sity function for the parameter space of the two models simul- 
taneously (since some of the free parameters are shared by both 
models). We used a Gaussian prior for the temperature of the 
star centered on the value measured by Shahbaz et al. (2022) 
during the disc state, T, = (6128 + 33) K. Gaussian priors have 
also been used for the parameters a, (centered on 0.2 + 0.2) and 
az (centered on —0.7 + 0.2), following the results of Bogdanov 
et al. (2018). We note that we also tried wider Gaussian priors for 
these parameters, resulting in wider posterior distributions with 
medians compatible with the results below within 1—2c. For all 
other parameters, we used uniform priors in very wide intervals 
to sample a large region of the parameter space (see last column 
of Table 8). This is of great importance to minimise the prob- 
ability that the actual source parameters fall outside the region 
sampled by the chosen priors. 


4.2.3. Results 


The results of the fit are reported in Table 8 and are summarized 
in the following paragraphs. All the parameters have been es- 
timated as the median of the marginal posterior distributions, 
with lo credible intervals coming from the 16—84th percentiles 
of the posterior distributions (see Fig. 10 for the histograms of 
the parameter samples). Most of the marginal posterior distri- 
butions of the parameters shown in Fig. 10 are nearly Gaussian, 
meaning that they are almost symmetrical. Fig. 10 also shows a 
strong correlation between Vsync and Fsync, and between Fo and 
a. This is however not surprising considering the scarce number 
of points at the lowest frequencies and the impossibility to con- 
strain quantities such as Vpreak and Vsync, due to numerous over- 
lapping components at infrared and optical frequencies. In par- 
ticular, we note that Fsync, Vsync and Qinin have a non-symmetrical 
posterior distributions, with long tails towards the upper or the 
lower limit of the chosen interval for the prior, depending on the 
parameter. However, the peak of their distributions is well de- 
termined (see Fig. 10). To estimate the goodness of the fit, we 


used the Bayesian p-value (reported in the last line of Table 8; 
Lucy 2016). The p-value is 0.78 for the high mode and 0.07 for 
the low mode, meaning that the results of the fit are indeed ac- 
ceptable (with a slight indication of overestimated and underesti- 
mated uncertainties for the high and low mode fit, respectively). 

The synchrotron emission caused by the shock has a higher 
flux density in the optical band than the average pulsed flux den- 
sity fraction of ~1% estimated by Papitto et al. (2019), which 
makes our scenario plausible. In particular, we find that the 
shock contributes ~3% of the emission in the optical band in 
the high mode, which is similar to the contribution of the opti- 
cally thin component of the compact jet emission. In the XMM- 
Newton/EPIC 0.3—10 keV range, the shock accounts for ~83% 
of the emission, while the jet only contributes ~17% of the flux 
in the high mode. 


4.3. Considerations on optical polarisation measurements 
and predictions of X-ray polarisation 


The average optical polarisation reported in this work is similar 
to levels previously measured for this source (Baglio et al. 2016; 
Hakala & Kajava 2018). Since the polarisation level is «196, it 
could have multiple origins. Given the lack of modulation of the 
Stokes parameters with the orbital period, it is unlikely that the 
source of the polarisation is scattering of radiation in any region 
of the system (such as the surface of the accretion disc), in con- 
trast to previous suggestions (Baglio et al. 2016). We identify 
two possible main emission mechanisms for the polarised emis- 
sion, both of which are non-thermal: optically thin synchrotron 
emission from the compact jet, or synchrotron emission at the 
shock between the pulsar wind and the accretion inflow, which 
also gives rise to the observed optical pulsations (as also sug- 
gested by Hakala & Kajava 2018). The jet flux density in the r- 
band is ~0.02 mJy, corresponding to a fractional contribution to 
the total flux in the same band of ~2.4% and ~2% in the high and 
low modes, respectively. In the high mode, where the polarisa- 
tion level measured in the R-band is 0.65+0.43%, this translates 
into a polarisation level of the jet emission of ~26% (or ~37% if 
the average polarisation level of 0.92% is considered instead). In 
the low mode, the jet emission is instead intrinsically polarised 
at a level of ~17%. Following Shahbaz (2019) and references 
therein, we estimate a level of ordering of the magnetic field lines 
of ~30% (considering the case of an optically thin synchrotron 
emission with a spectral slope of ~-0.78, as obtained from our 
fit and averaging across modes; see Table 8). This indeed indi- 
cates a highly ordered magnetic field and suggests that tMSPs 
might have a more ordered magnetic field at the base of the jet 
than other NS and BH X-ray binaries (XRBs) (Russell & Fender 
2008; Baglio et al. 2014; Russell et al. 2016; Baglio et al. 2020). 
Based on our picture, this could be ascribed to the presence of 
a continuous, steady flow of plasma near the pulsar, which may 
help to maintain a more ordered magnetic field structure at the 
base of the jet. 

Assuming that the observed optical polarisation is due to this 
emission component and that the same population of accelerated 
electrons in the jet produces optically thin synchrotron radiation 
at all frequencies above the break frequency, we can predict a po- 
larisation level in the X-ray band (in the 2-8 keV energy range of 
the Imaging X-ray Polarimetry Explorer, /XPE; Weisskopf et al. 
2016; Soffitta et al. 2021) due to the jet of «17406 in the low mode, 
where the jet contributes to ~100% of the flux according to our 
SED modeling. 

Alternatively, the observed optical polarisation might origi- 
nate from synchrotron radiation at the region of shock between 
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Table 8. Results of the model fit of the SED in low and high modes. The median values of the posterior distributions of the model parameters are 
reported, along with lower and upper uncertainties. These uncertainties correspond to the 15.9th and 84.1th percentiles of the posterior distributions 
for each parameter. The last column reports the prior bounds that were used for all parameters; AN (u, o`) is a normal distribution centered in u and 


with variance o? 


, and (a, b) is a uniform distribution in the (a, b) interval. 


Component Parameter Posterior Prior bounds 
Low mode High mode 
Irradiated star T, (K) 6141.05 72271 ~ N(6128,33)* 
Accretion disc — Log [Finopt/uv / 1 cm] 9.49 + 0.02 Logrinopyuv ~ U(7, 12) 
Compact jet Log [break / 1 Hz] 13.40 + 0.26 Log Vbreak ~ U(10, 14) 
Q1 0.18 + 0.06 ~ N(0.2, 0.2) 
a —0.78 + 0.04 ay ~ N(-0.7,0.2)' 
Log [Fo / 1 mJy] —0.59*0.5 Log Fo ~ U(-2, 0) 
Shock emission Log [Vsync / 1 Hz] — 14. 3404 Log Vsyne ~ U(12, 16) 
Ü'sync — —0.61: 0.01  @syne ~ U(—0.9, —0.4) 
Log [F'sync /1 mJy] = =]; 1140 2 Log Fsync = Uu(-3, 2) 
Discrete ejecta ^ ain —0. 19 2 - Qhin ~ U(-1.5, 0) 
Log [Fininio/ 1 mJy] —1.03 + 0.05 —  LogFmino ~ U(-S, 2) 
Bayesian p-value 0.07 0.78 


Notes. * 


the pulsar wind and the accretion flow. In this scenario, polar- 
isation should only be detected during the high mode. Unfor- 
tunately, the significance of our polarisation results obtained in 
the high and low modes separately is not high enough to ver- 
ify this, as the results are compatible with each other within 1c 
(0. 36*? 270b in the low mode and 0.65+0.43% in the high mode). 
The fractional flux contribution of the shock in the r-band ac- 
cording to the SED modeling in the high mode is ~4.8%, which 
translates into a linear polarisation of the emission of the shock 
of ~15-20%, with peaks of up to 65-90% at the pulsations (in 
fact, the fraction of pulsed radiation is typically of ~1%; Pa- 
pitto et al. 2019). This prediction could be tested using next- 
generation optical high-time resolution polarimetric instruments 
(Ghedina et al. 2022). In the X-rays instead, we predict a frac- 
tional flux contribution of the shock emission of ~90% in the 
IXPE energy range, which converts into an intrinsic linear po- 
larisation level of «12-1796. This could be detected in an IXPE 
observation of ~460 ks (assuming J1023 spends 70% of the time 
in the high mode). 


In both scenarios, the validity of these predictions depends 
on the assumption that the optical polarisation is solely caused 
by synchrotron radiation from the jet or shock. However, if the 
observed polarisation is due to a combination of these factors or 
if other processes are at play that are not being considered (such 
as scattering on the disc surface contributing to the modulation 
with a low SNR in our measurements), these predictions may 
need to be revised. 


4.4. The millimetre continuum emission of J1023 and other 
X-ray binary systems 


We have reported the detection of millimetre continuum emis- 
sion for the first time from J1023 (and from a tMSP in gen- 
eral) using ALMA observations at «100 GHz. The millimetre 
emission is unresolved, implying a size for the emission region 
smaller than ~ 6 x 10/5 cm (this corresponds to ~ 7 x 104 
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Value from Shahbaz et al. (2022). * Value from Bogdanov et al. (2018). 


times the orbital separation of the system; Archibald et al. 2013). 
Millimetre/sub-millimetre continuum emission has been firmly 
detected from a number of transiently accreting stellar-mass 
black holes (Paredes et al. 2000; Ogley et al. 2000; Fender et al. 
2000, 2001; Ueda et al. 2002; Russell et al. 2013b; van der Horst 
et al. 2013; Tetarenko et al. 2015, 2017, 2019; Russell et al. 
2020; de Haas et al. 2021; Koljonen & Hovatta 2021; Tetarenko 
et al. 2021; Díaz Trigo et al. 2021) and in three accreting NSs so 
far: the persistent systems 4U 1728—34 and 4U 1820—30 (Díaz 
Trigo et al. 2017) and the transient accreting MSP Aql X-1 (Diaz 
Trigo et al. 2018). In these cases, however, such emission has 
been detected during accretion states that are seemingly differ- 
ent from the X-ray sub-luminous disc state attained by J1023 at 
the epoch of our ALMA observations, as detailed below. 

In black-hole systems, millimetre continuum emission is typ- 
ically detected during the rising phase of an X-ray outburst 
or along the outburst decay, during the "hard" X-ray spectral 
state. In this state, the X-ray emission is typically dominated by 
a power-law component that is believed to arise from Comp- 
ton up-scattering of soft disc photons by a population of hot 
thermal electrons located either in the innermost region of the 
accretion flow, or in an extended cloud above the disc (for a 
recent review, see Motta et al. 2021). At low frequencies, a 
flat to slightly inverted optically thick spectrum (a > 0) has 
been observed from radio through sub-millimetre frequencies 
and above and has been associated with partially self-absorbed 
synchrotron radiation from electron populations in different re- 
gions along a compact jet. This spectrum becomes optically thin 
(a < 0) typically around infrared frequencies, resulting in a 
spectral break (Russell et al. 2013a,c). On the other hand, at 
very high mass accretion rates (typically close to the outburst 
peak), these systems can be in a “soft” X-ray spectral state where 
the X-ray emission is dominated by a thermal component pro- 
duced in the hot inner regions of the accretion disc (see Motta 
et al. 2021 and references therein). In the state that immedi- 
ately precedes or follows the soft state during the outburst (the 
so-called 'soft-intermediate state"), the optically-thick jet emis- 
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Fig. 10. Corner plot showing the posterior probability distributions of the parameters obtained from the MCMC sampling algorithm. The solid 
lines in the panels along the diagonal represent the probability density of each parameter. The vertical dashed lines indicate the 16th, 50th, and 
84th percentiles of the distributions. The median of the parameter posterior and the 1c error bars are displayed above each panel. The off-diagonal 
panels show the two-dimensional posterior distributions with contour lines indicating the 1-, 2-, and 3-c equivalent bounds. 


sion is observed to be quenched altogether, while discrete jet 
ejecta characterized by an optically thin steep spectrum (œ < 0) 
in the radio band can be launched (Fender et al. 2004; Rus- 
sell et al. 2020) and be accompanied by strong radio variabil- 
ity and bright multi-frequency flaring activity (even in the sub- 
millimetre band; Tetarenko et al. 2019). 

The picture for NS systems is more diverse (see van den Ei- 
jnden et al. 2021 for a recent overview). All three of the above- 
mentioned NS systems were observed in the millimetre band 
with ALMA, along with observations at other frequencies, in 
regimes of high mass accretion rates (Ly > 5 x 10? erg s^!). 
4U 1728—34 was observed during a transition from a hard to 
a soft state. The low-frequency spectrum was compatible with 
synchrotron emission from a compact jet and also revealed the 


presence of a break from optically thick to optically thin syn- 
chrotron radiation at frequencies in the range 10? — 10!^ Hz. 
In the case of 4U 1820—30, enhanced millimetre emission was 
detected during a soft X-ray state (the first time millimetre emis- 
sion has been observed for an accreting NS in the soft state). The 
spectrum from the radio to the sub-millimetre bands was slightly 
inverted and showed no evidence for the presence of a spectral 
break (Díaz Trigo et al. 2017). This indicates emission from a 
compact jet and supports earlier findings that jets in NS systems 
do not all appear to be completely suppressed in soft accretion 
states (Migliari et al. 2004). Aql X-1 was observed five times 
using ALMA over a time span of ~1.5 months during the decay 
phase of an X-ray outburst in 2016, covering multiple transi- 
tions between distinct accretion states. The spectrum from the 
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Fig. 11. Broadband (radio to X-rays) SEDs of J1023 during the high and low modes, along with the best-fitting models. Both SEDs were corrected 
for extinction effects. The red-dashed, blue-dotted, purple-dashed-dotted and orange solid lines represent the contributions of the accretion disc, 
companion star, discrete ejection and shock region, respectively. The green dashed line represents the compact jet, whose contribution does not 


vary between the two modes. The green solid line is the sum of all contributions in both SEDs. The ratio between the data points and the best-fitting 
model is shown at the bottom of each panel. 


Table 9. Millimetre and X-ray luminosities of accreting compact objects in X-ray binaries. Values are reported for those systems that have been 
clearly detected in the millimetre/sub-millimetre bands, have a known distance, and have simultaneous or nearly simultaneous millimetre/sub- 
millimetre and X-ray observations. Systems above the horizontal line are black-hole systems, while those below the line are NS systems. 


Source Distance  Millimetre luminosity? X-ray luminosity? Reference 
(kpc) (erg/s) (erg/s) 
MAXIJ1659-152 6 (7x1)x10?! (1.389+0.001)x1037 van der Horst et al. (2013) 
MAXIJ1659-152 6 (5x1)x10?! (4.60x0.01)x10?7 van der Horst et al. (2013) 
MAXIJ1659-152 6 «2.6x10?! (4.8120.01))x10?7 van der Horst et al. (2013) 
MAXIJ1659-152 6 «21x10?! (4.67x0.02)x10?7 van der Horst et al. (2013) 
V404 Cygni 2.39 (3.64+0.06)x10°° (8+2)x10%4 Tetarenko et al. (2019) 
V404 Cygni 2.39 (9.42-0.6)x 1079 (1.5+0.1)x 1034 Tetarenko et al. (2019) 
MAXIJ1535-571 4.1 (4.70.2)x10?? (2.00+0.01)x10°8 Russell et al. (2020) 
MAXIJ1535-571 4.1 (1.0+0.2)x 103? (8.6+0.1)x1038 Russell et al. (2020) 
GX 339-4 8 (1.53+0.08)x10?! (1.3540.08)x10?6 de Haas et al. (2021) 
MAXIJ1820+070 2.96 (9.651+0.004)x 103! (3.50+0.01)x10°7 Tetarenko et al. (2021) 
GRS 1915+105 8.6 22x10?! (2.07+0.03)x10°7 Koljonen & Hovatta (2021) 
GRS 1915+105 8.6 2.6x10?! (2.69:0.01)x10?7 Koljonen & Hovatta (2021) 
XTE J1118+480 1.7 6.9x109?9 2.4x10%5 Fender et al. (2001) 
4U 1820-30* 6.4 (1.67+0.09)x10™ (3.22+0.03)x10°7 Diaz Trigo et al. (2017) 
Aql X-1 5.2 (2.06+0.08)x 103? (1.73+0.01)x10?7 Díaz Trigo et al. (2018) 
Aql X-1 5.2 (1.31+0.08)x10°° (9.2+0.1)x10°” Diaz Trigo et al. (2018) 
Aql X-1 52 (2.44+0.09)x10°° (6.65+0.07)x10°° Diaz Trigo et al. (2018) 
PSR J1023+0038 (mean) 1.37 (3.2+0.2)x 1078 (1.7+0.2)x 108 This work 
PSR J1023+0038 (high) 1.37 (2.0+0.3)x 1078 (2.40+0.04)x 1033 This work 
PSR J1023+0038 (low) 1.37 (3.020.5)x1028 (3.620.4)x10?? This work 


Notes. ? Millimetre luminosities are evaluated at a frequency of 100 GHz. A flat spectrum was assumed to estimate the luminosity 
for those cases where the flux density is available at other frequencies in the literature, except for MAXI J1820--070 (@=0.25) and 
4U 1820-30 (a0.13). 
^ X-ray luminosities are evaluated over the 1-10 keV energy range and are derived from Swift/XRT observations except for the 
case of J1023, for which values are taken from Table 5. Spectra were extracted using the Swifi/XRT data products generator (Evans 
et al. 2009) and fitted using an absorbed power-law model. For the second pair of simultaneous observations of GRS 1915+105, 
the X-ray luminosity is evaluated using NJCER obs ID 2596010301. 

* 'The millimetre and X-ray observations were separated by 4 days. 
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Fig. 12. Millimetre and X-ray luminosities for accreting compact objects in X-ray binaries. BH systems are marked in black. The adopted values 
for the luminosities are listed in Table 9. The black dashed line and the gray shaded area show the best-fitting slope for the BH correlation and the 
region on the plane at less than 3° from the best-fitting slope, respectively (see text for details). 


radio to the millimetre bands was again compatible with emis- 
sion from a compact jet and, for the first time in an accreting NS 
in a low-mass X-ray binary, a break in the SED was detected at 
a frequency that varies dramatically depending on the accretion 
state. Specifically, the break frequency initially decreased from 
~100 GHz to less than ~5 GHz during a transition from a hard 
to a soft state, and then increased again up to a frequency within 
the range ~30-100 GHz during a transition to a hard state at later 
stages of the outburst decay (Díaz Trigo et al. 2018). 


Unlike all the cases mentioned above, the ALMA observa- 
tions of J1023 were performed at an epoch when the system did 
not display any changes in its spectral state and was lingering 
in a sub-luminous X-ray state («10?? erg s~!). It is worth noting 
that both the average and the mode-resolved X-ray spectra of 
J1023 have always been hard and have not shown any substan- 
tial changes over the past ten years or so during the sub-luminous 
X-ray state. 


We can compare the millimetre and X-ray luminosities of 
J1023 with those measured for black-hole and NS systems using 
quasi-simultaneous millimetre and X-ray observations. These lu- 
minosities are reported in Table9 and plotted in Fig. 12. De- 
spite the relatively small number of existing millimetre/sub- 
millimetre data for X-ray binaries, Fig. 12 shows that millime- 


tre continuum emission has been detected in systems that at- 
tain a wide range of X-ray luminosities (-107?—10?? erg s^!) and 
that J1023 currently holds the record as the faintest X-ray binary 
firmly detected in the millimetre band. 

The mechanism responsible for the X-ray emission in the 
two modes in J1023 is clearly different from that at work in 
black-hole systems in the hard state, so a proper comparison 
cannot be made. With this in mind, we assessed how the po- 
sition of J1023 on the Lmm — Lx plane compares with those of 
black-hole systems by parameterizing the dependence of Linm 
upon Z4 as a relation of the form LogLmm — LogLmmave = 
a + B(LogLx — LogLx ay.) ^. We then investigated the shape of 
the Lim — Lx correlation for black-hole systems by performing 
a Bayesian-based MCMC sampling of the black-hole data (in- 
cluding upper limits) using the linear regression algorithm Lm- 
MIX ERR (Kelly 2007). There are three free parameters: the in- 
tercept a, the slope B and oo, a parameter that accounts for an 
intrinsic random (Gaussian) scatter of the values around the best- 
fit model. We calculated the median values of the fit parame- 
ters from 10^ draws from the posterior distribution, and deter- 
mined the lo credible intervals of these parameters using the 


E Lunae = 103193 erg s^! and Lxave = 1038.03 erg s^! are the average 


values for the detections in our BH sample. 
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16-84th percentiles of the posterior distributions. We derived 
a = —0.04 + 0.15, B = 0.46 + 0.09 and op = 0.18 + 0.13. The 
best-fitting slope is consistent with that derived for the Lradio — Lx 
correlation for black-hole systems (Gallo et al. 2018) within the 
uncertainties. Fig. 12 shows the best-fitting slope (black dashed 
line) and the region on the plane where the values of Lj and 
Lx are located at a distance of less than 30 from the best-fitting 
slope (gray shaded area). This figure shows that, during the low 
mode, J1023 follows the same correlation as black-hole systems 
in the hard state (Bogdanov et al. 2018). 


5. Summary and Conclusions 


We have presented the results of the largest multi-wavelength 
campaign ever conducted on the tMSP PSR J1023--0038. This 
campaign took place in June 2021, while the system was in an 
active low-luminosity X-ray disc state. Our key findings can be 
summarized as follows: 

e During the first night of observations, flaring activity was 
detected at X-ray and optical wavelengths, followed by mode- 
switching behavior in X-rays. The UV emission was signifi- 
cantly fainter in low mode, while the radio continuum emission 
was brighter. No evidence of pulsed radio emission was detected 
in either mode. Our observations provide the most stringent lim- 
its on the pulsed radio flux density to date. 

e During the second night observations, a low-significance 
variation in the optical polarised signal was detected during 
mode changes, with a slightly higher polarisation level in the 
high mode compared to the low mode. The millimetre emis- 
sion did not strictly follow the X-ray low modes, except for two 
brief flares. One of these flares occurred at the high-to-low mode 
switch, while the exact onset time of the other flare could not be 
determined due to a lack of observational coverage. 

e These results allow us to draw a physical picture of high- 
low mode switches in J1023, involving a rotation-powered pul- 
sar, an accretion disc, and discrete mass ejections on top of a 
relativistic compact jet. During the high mode, synchrotron ra- 
diation at the shock front between the pulsar wind and the inner 
accretion flow produces most X-ray emission, as well as X-ray, 
UV, and optical pulsations. The switch to low mode is triggered 
by discrete mass ejections, resulting in a drop of the X-ray flux 
and quenching of pulsations. During the low mode, the pulsar 
wind still penetrates the accretion disc and launches a compact 
jet, contributing to the multiband emission in the same way as in 
the high mode. Eventually, the flow from the disc refills the inner 
regions outside NS light cylinder, leading to increased emission 
and pulsations at X-ray, UV, and optical frequencies. The sys- 
tem switches back to the high mode. This scenario is supported 
by the modelling of the SEDs extracted separately during the 
high and low modes. 

e Relativistic pulsar winds can shape the structure and dy- 
namics of accretion flows in X-ray binary systems containing 
quickly spinning NSs, and can be used to probe emission mech- 
anisms in tMSPs and other X-ray binary systems. 

e tMSPs and black-hole systems share intriguing similari- 
ties in their phenomenological properties, emphasizing the need 
for further research to deepen our understanding of accretion 
physics in compact objects. 
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Appendix A: Searches for counterparts to 
ALMA J102349.1--003824 


We searched for a counterpart to ALMA J102349.1+003824 at 
other wavelengths using both simultaneous and archival ob- 
servations of the field. No UV emission was detected at the 
source position in the Swift/UVOT observations simultaneous 
with ALMA, either in the images collected separately during 
the three snapshots, or in the stack of these images. We set a 
30 upper limit on the magnitude of a putative counterpart of 
UVM2>21 from the stacked image (total exposure of ~2.4 ks). 
A few weeks before the ALMA observations were acquired, the 
field of view of ALMA J102349.1--003824 was also observed 
in the optical (g', r’ and i’ SDSS filters) with the Las Cumbres 
Observatory (LCO) network of telescopes. No counterpart was 
detected in these data. We set the following 30 upper limits on 
the optical magnitude using 60-s images acquired with the 1m 
telescope of the LCO network located in South Africa on May 
8, 2021: i’ > 17.8, r' > 17.9, g' > 18.6. There is also no opti- 
cal/NIR counterpart in either the third data release (DR3) from 
the Gaia space telescope or in the 2MASS point source catalog. 
ALMA J102349.1+003824 was also not detected in the VLA ob- 
servation carried out a few weeks earlier, down to a 3-c- flux 
density upper limit of 12 uJy at 10 GHz. 

Several observations of the field have been performed 
in the soft X-ray band over the past years, mainly us- 
ing Swift and XMM-Newton. However, given the X-ray lu- 
minosity of J1023 and its small angular distance from 
ALMA J102349.1+003824, the point-spread function of the X- 
ray instruments on board these observatories encompasses the 
position of ALMA J102349.1+003824, preventing a conclusive 
assessment of the presence or absence of an X-ray counterpart 
to this source. We then analyzed a 83-ks long Chandra/ACIS-S 
observation performed in March 2011 (Obs. ID 11075) to search 
for X-ray emission at the position of ALMA J102349.1+003824. 
We reprocessed the data using the task CHANDRA_REPRO Within the 
ciao software and ran the wavpErECT source detection tool. No 
evidence for an X-ray counterpart to ALMA J102349.1--003824 
was found. Using sRcrLUx, we set a 3c upper limit on the net 
count rate of 1.3 x 107^ counts s^! over the 0.5-7.0 keV energy 
range. Assuming an absorbed power-law spectrum with an ab- 
sorption column density of Ny = 5 x 10? cm"? (i.e., the one 
expected within the Galaxy in the direction of the source: ?) and 
photon index in the range I = 1 — 4, the limit above translates 
into an unabsorbed flux of < (2-4)x10- P? erg cm"? s! (depend- 
ing on the assumed power-law slope) in the 0.5—7 keV range. 
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